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Abstract 

Motility initiation in crawling cells requires transformation of a symmetric state into a polarized state. In contrast, 
motility arrest is associated with re-symmetrization of the internal configuration of a cell. Experiments on kera- 
tocytes suggest that polarization is triggered by the increased contractility of motor proteins but the conditions of 
re-symmetrization remain unknown. In this paper we show that if adhesion with the extra-cellular substrate is suffi¬ 
ciently low, the progressive intensification of motor-induced contraction may be responsible for both transitions: from 
static (symmetric) to motile (polarized) at a lower contractility threshold and from motile (polarized) back to static 
(symmetric) at a higher contractility threshold. Our model of lamellipodial cell motility is based on a ID projection of 
the complex intra-cellular dynamics on the direction of locomotion. In the interest of analytical transparency we also 
neglect active protrusion and view adhesion as passive. Despite the unavoidable oversimplifications associated with 
these assumptions, the model reproduces quantitatively the motility initiation pattern in fish keratocytes and reveals a 
crucial role played in cell motility by the nonlocal feedback between the mechanics and the transport of active agents. 
A prediction of the model that a crawling cell can stop and re-symmetrize when contractility increases sufficiently far 
beyond the motility initiation threshold still awaits experimental verification. 


1. Introduction 

The ability of cells to self-propel is essential for many biological processes. In the early life of an embryo, stem 
cells move to form tissues and organs. During the immune response, leukocytes migrate through capillaries to attack 
infections. Wound healing requires the motion of epithelial cells. While the biochemistry of such motility is rather 
well understood, the underlying mechanics of active continuum media is still in the stage of development (Bray, 2000; 
Mogilner, 2009; Carlsson and Sept, 2008; Joanny and Frost, 2011; Adler and Givli, 2013; Ziebert and Aranson, 2013; 
Giomi and DeSimone, 2014; Recho et al., 2014; Cox and Smith, 2014). 

At a very general level, a cell can be viewed as an elastic ‘bag’ whose interior is separated from the exterior by a 
bi-layer lipid membrane. The membrane is attached from inside to a thin cortex - an active muscle-type layer main¬ 
taining the cell’s shape. The interior is filled with a passive medium, the cytosol, where all essential cell organelles 
are immersed. The active machinery inside the cytosol ensuring self-propulsion is contained in the cytoskeleton: a 
perpetually renewed network of actin filaments cross-linked by myosin motors that can inflict contractile stresses. The 
cytoskeleton can be mechanically linked to the cell exterior through adhesion proteins (Alberts et al., 2002). 

The elementary mechanisms responsible for the steady crawling of keratocytes (flattened cells with fibroblastic 
functions) have been identified (Abercrombie, 1980; Bell, 1984; Stossel, 1993; Bellairs, 2000). The advance starts 
with protrusion driven by active polymerization of the actin network in the frontal area of the cell (the lamellipodium) 
with a simultaneous formation of adhesion clusters at the advancing edge. After the adhesion of the protruding part 
of the cell is secured, the cytoskeleton contracts thanks to the activity of myosin motors. This contraction leads to the 
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detachment at the rear and lessening of the actin network through de-polymerization. All these active phenomena are 
driven by ATP hydrolysis and are highly synchronized which allows the cell to move with a stable shape and relatively 
constant velocity (Barnhart et al., 201 1). 

The initiation of such motility requires a polarization of the cell which is a process that discriminates the leading 
edge from the trailing edge. The implied symmetry breaking turns a symmetric stationary configuration of a cell 
into a polar motile configuration. While both contraction and protrusion contribute to steady state cell migration, 
contraction appears to be the dominating mechanism of polarization: it has been shown experimentally that motility 
initiation in keratocytes may be triggered by raising the contractility of myosin (Verkhovsky et al., 1999; Csucs et al., 
2007; Lombardi et al., 2007; Yam et al., 2007; Vicente-Manzanares et al., 2009; Poincloux et al., 2011). It is also 
known that cells may self propel by contraction only (Keller et al., 2002). 

In physical terms, the contraction-driven polarization/motility is performed by ‘pullers’ (contractile agents) while 
‘pushers’ (protrusive agents) remain largely disabled. Some numerical models suggest that the relative role of ‘push¬ 
ers’ and ‘pullers’ in cellular motility may be tightly linked to the task to be performed (Simha and Ramaswamy, 2002; 
Saintillan and Shelley, 2012) and even to the nature of the cargo (Recho and Truskinovsky, 2013). However, it is still 
not fully understood why in case of keratocytes the motility initiation is primarily contraction-driven. In contrast to 
motility initiation, the reciprocal process of motility arrest is associated with re-symmetrization and such symmetry 
recovery is a typical precursor of cell division (Stewart et al., 2011 ; Lancaster et al., 2013; Lancaster and Baum, 2014). 
It is not yet clear whether re-symmetrization is also predominantly contraction-driven and if yes, whether it requires 
contractility reduction or contractility increase beyond the motility initiation threshold. It is, however, known that 
some cells can switch from static to motile state as a result of a decrease in the level of contractility (Liu et al., 2010; 
Hur etal., 2011). 

A large variety of modeling approaches targeting cell motility mechanisms can be found in the literature, see 
the reviews by Rafelski and Theriot (2004); Carlsson and Sept (2008); Mogilner (2009); Wang et al. (2012). In some 
models, the actin network is viewed as a highly viscous dLCtive fluid moving through a cytoplasm by generating internal 
contractile stresses (Alt and Dembo, 1999; Oliver et al., 2005; Herant and Dembo, 2010; Kimpton et al., 2014). In 
other models, the cytoskeleton is represented by an active gel whose polar nature is modeled in the framework of 
the theory of liquid crystals (Kruse et al., 2005; Joanny et al., 2007; Julicher et al., 2007; Joanny and Prost, 2011; 
Callan-Jones and Julicher, 2011). The active gel theory approach, which we basically follow in this study without 
an explicit reference to local orientational order, was particularly successful in reproducing rings, asters, vortices and 
some other sub cellular structures observed in vivo (Doubrovinski and Kruse, 2007; Sankararaman and Ramaswamy, 
2009; Doubrovinski and Kruse, 2010; Du et al., 2012). At sufficiently fast time scales, the cytoskeleton can be also 
modeled as an active solid with a highly nonlinear scale-free rheology (Broedersz and MacKintosh, 2014; Pritchard 
et al., 2014). 

Various specific sub-elements of the motility mechanism have been subjected to careful mechanical study. Thus, it 
was shown that in some cases the plasmic membrane with an attached cortex can be viewed as a passive elastic surface 
and modeled by phase field methods allowing one to go smoothly through topological transitions (Wang et al., 2012; 
Giomi and DeSimone, 2014). In other cases, the membrane may also play an active role, for instance, an asymmetric 
distribution of channels on the surface of the membrane can be responsible for a particular mechanism of cell motility 
relying on variation of osmotic pressure (Stroka et al., 2014). While most models assume that the cell membrane 
interacts with the exterior of the cell through passive viscous forces, active dynamics of adhesion complexes has 
recently become an area of intense research driven in part by the finding of a complex dependence of the crawling 
velocity on the adhesive properties of the environment (DiMilla et al., 1991; Novak et al., 2004; Deshpande et al., 
2008; Gao et al., 2011; Lin et al., 2008; Ronan et al., 2014; Lin, 2010; Ziebert and Aranson, 2013). The account 
of other relevant factors, including realistic geometry, G-actin transport, Rac/Rho-regulation, etc., have led to the 
development of rather comprehensive models that can already serve as powerful predictive tools (Rubinstein et al., 
2009; Wolgemuth et al., 2011; Tjhung et al., 2012; Giomi and DeSimone, 2014; Barnhart et al., 2015). 

The more focused problem of finding the detailed mechanism of motility initiation, is most commonly ad¬ 
dressed in the framework of theories emphasizing polymerization-driven protrusion (Mogilner and Edelstein-Keshet, 
2002; Dawes et al., 2006; Bernheim-Groswasser et al., 2005; Schreiber et al., 2010; Campas et al., 2012; Hodge 
and Papadopoulos, 2012). With such emphasis on ‘pushers’, spontaneous polarization was studied by Kozlov and 
Mogilner (2007); Callan-Jones et al. (2008); John et al. (2008); Hawkins et al. (2009); Hawkins and Voituriez (2010); 
Doubrovinski and Kruse (2011); Blanch-Mercader and Casademunt (2013). In Banerjee and Marchetti (201 1); Ziebert 
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et al. (2012) and Ziebert and Aranson (2013), polarization was interpreted as a result of an inhomogeneity of adhesive 
interactions. Yet another group of authors have successfully argued that cell polarity may be induced by a Turing-type 
instability (Mori et ah, 2008; Altschuler et ah, 2008; Vanderlei et ah, 2011; Jilkine and Edelstein-Keshet, 2011). Such 
a variety of modeling approaches is a manifestation of the fact that very different mechanisms of motility initiation 
are engaged in cells of different types. 

The observation that contraction may be the leading factor behind the polarization of keratocytes has been broadly 
discussed in the literature. It was realized that active contraction creates an asymmetry-amplifying positive feedback 
because it causes actin flow which in turn carries the regulators of contraction (Kruse et al., 2003; Ahmadi et al., 2006; 
Salbreux et al., 2009; Recho et al., 2013; Barnhart et al., 2015). In constrained conditions such positive feedback gen¬ 
erates peaks in the concentration of stress activators (myosin motors) (Bois et al., 2011; Howard et al., 2011) and this 
patterning mechanism was used to model polarization induced by angular cortex flow (Hawkins et al., 2009, 2011). 
Closely related heuristic models of the Keller-Segel type (Perthame, 2008) describing symmetry breaking and local¬ 
ization were independently proposed by Kruse and Julicher (2003); Calvez et al. (2010). In all these models, however, 
the effect of contraction (pullers) was obscured by the account of other mechanisms, in particular, polymerization 
induced protrusion (pushers), and the focus was on generation of internal flow and the resulting pattern formation, 
rather than on the problem of ensuring steady translocation of a cell. 

This shortcoming was overcome in more recent models of contraction-induced polarization relying on splay insta¬ 
bility in an active gel (Tjhung et al., 2012; Giomi and DeSimone, 2014; Tjhung et al., 2015). In these model, however, 
‘pushers’ were not the only players, in particular, polarization was induced by a local phase transition from non-polar 
to polar gel. In Callan-Jones and Voituriez (2013), the motility initiation was attributed to a contraction-induced in¬ 
stability in a poro-elastic active gel permeated by a solvent. Here again the non-contractile active mechanism was 
involved as well and therefore the domineering role of contraction could not be made explicit. 

The goal of the present paper is to focus on the special role of bare contraction in symmetry breaking processes 
by studying a minimalistic, analytically transparent model of motility initiation in a segment of an active gel. Fol¬ 
lowing previous work, we exploit the Keller-Segel mechanism, but now in difree boundary setting, and show that the 
underlying symmetry breaking instability is fundamentally similar to an uphill diffusion of the Cahn-Hilliard type. In 
contrast to most previous studies, our contraction driven translocation of a cell is caused exclusively by the internal 
flow generated by molecular motors (pullers) and no other active agents are involved. Each ‘puller’ contributes to the 
stress field and simultaneously undergoes biased random motion resulting in an uphill diffusion along the correspond¬ 
ing stress gradient. In other words our ‘pullers’ (active cross-linkers) use the continuum environment (passive actin 
network) as a medium through which they interact and self-organize. 

We emphasize that the contraction mechanism of polarization and motility (Recho et al., 2013, 2014) is concep¬ 
tually very close to chemotaxis, however, instead of chemical gradients, the localization and motility is ultimately 
driven by the self-induced mechanical gradients. More specifically, the pullers propel the passive medium by inflict¬ 
ing contraction which creates an autocatalytic effect since the pullers are themselves advected by this medium (Mayer 
et al., 2010). The inevitable build up of mechanical gradients in these conditions is limited by diffusion which resists 
the runaway and provides the negative feedback. After the symmetry of the static configuration is broken in the con¬ 
ditions where matter can circulate, the resultant contraction-driven flow ensures the perpetual renewal of the network 
and then frictional interaction with the environment allows for the steady translocation of the cell body. 

The next natural question is how such steady translocations can be halted. For instance, if motility initiation is 
contraction-driven, can motility arrest be also contraction driven and what a steadily moving cell can do in order to 
stop and symmetrize? Several computational models provided an indication that motility initiation and motility arrest 
may be related to a re-entrant behavior of the same branch of motile regimes (Kruse and Julicher, 2003; Tjhung et al., 
2012; Recho et al., 2013; Giomi and DeSimone, 2014). To make the link between motility initiation and motility 
arrest more transparent we study in this paper an analytically tractable problem which captures the complexity of the 
underlying physical phenomena. While most of the elements of the proposed model (Recho et al., 2013, 2014) have 
been anticipated by some comprehensive computational approaches (e.g. Rubinstein et al., 2009), it was previously 
not apparent that the initiation of motility, steady translocation and the arrest of motility can be all captured in such a 
minimal setting. 

Our model of lamellipodial cell motility is based on a ID projection of the complex intra-cellular dynamics 
on the direction of locomotion. In the interest of analytical transparency, we decouple the dynamics of actin and 
myosin by assuming infinite compressibility of the cross-linked actin network (Julicher et al., 2007; Rubinstein et al.. 
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2009). To ensure that the crawling cell maintains its size, we introduce a simplified cortex/osmolarity mediated 
quasi-elastic interaction between the front and the back of the self-propelling fragment (Banerjee and Marchetti, 
2012; Barnhart et al., 2010; Du et al., 2012; Loosley and Tang, 2012); a comparison of such mean field elasticity 
with more conventional bulk elasticity models can be found in Recho and Truskinovsky (2013). We remark that the 
coupling between the front and the rear of the fragment may also have an active origin resulting from different rates of 
polymerization and depolymerization at the extremities of the lamellipodium (Recho and Truskinovsky, 2013; Etienne 
et al., 2014). In other respects we neglect active protrusion (pushers). We also view adhesion as fully passive. 

Despite the unavoidable oversimplifications associated with these assumptions, we show that our model repro¬ 
duces quantitatively the motility initiation pattern in fish keratocytes and reveals a crucial role played in cell motility 
by the nonlocal feedback between the mechanics and the transport of active agents. It also provides compelling evi¬ 
dence that both, the initiation of motility and its arrest, may be fully controlled by the average contractility of motor 
proteins. 

More precisely, we show that the increase of contractility beyond a well defined threshold leads to a bifurcation 
from a static symmetric solution of the governing system of equations (of Keller-Segel type) to an asymmetric trav¬ 
eling wave (TW) solution corresponding to steadily moving cells. While several TW regimes may be available at 
the same value of parameters, we show that stable TW solutions localize motors at the trailing edge of the cell in 
agreement with observations (Verkhovsky et al., 1999; Csucs et al., 2007; Lombardi et al., 2007; Yam et al., 2007; 
Vicente-Manzanares et al., 2009; Poincloux et al., 2011). Moreover, we show that if adhesion with the extra-cellular 
substrate is sufficiently low, the increase of motor-induced contraction may induce transition from the steady state 
TW solution back to a static solution. This re-symmetrization transition, leading to the motility arrest, can be directly 
associated with the behavior of keratocytes prior to cell division and our model shows that such re-entrant behavior 
can be ensured by ‘pullers’ without any engagement of either active protrusion or liquid crystal elasticity. 

The paper is organized as follows. In Section 2, we present a discrete “model of a model” which conveys the main 
ideas of our approach in the simplest form. In Section 3, we develop a continuum analogue of the discrete model, 
study its mathematical structure and pose the problem of finding the whole set of TW solutions incorporating both 
static and motile regimes. In Section 4, all static solutions of the TW problem are found analytically. In Section 5, we 
study the fine structure of multiple bifurcations producing motile solutions from the static ones and identify parametric 
regimes when these bifurcations become re-entrant. In Section 6, we investigate numerically the initial value problem 
which allows us to qualify some of the motile TW solutions as attractors. The reconstruction the background turnover 
of actin, which takes place in our model without active protrusion at the leading edge, is discussed in Section 7. 
In Section 8, we demonstrate that our model can quantitatively match the experiments carried on keratocytes. The 
last Section highlights our main conclusions and mentions some of the unsolved problems; three appendices contain 
material of technical nature. 

Some of the results of this paper have been previously announced in two pre-publications (Recho et al., 2013, 
2014) but without any details. In addition to providing a necessary background for Recho et al. (2013, 2014), here 
we develop a new discrete model, investigate the nonlocal nature of the coupling between mechanics and diffusion 
of active agents, give a thorough analysis of the static regimes, study the bifurcation points by using the Lyapunov- 
Schmidt reduction technique, investigate the non-steady problem numerically, generalize the model to account for 
nonlinear dependence of contractile stress on motor concentration and provide a detailed quantitative comparison of 
the model with experiment. 

2. The discrete model 

Our point of departure is a conceptual discrete model elucidating the mechanism of contraction-driven crawling 
and emphasizing the role of symmetry breaking in achieving the state of steady self propulsion. This ’’model of a 
model” allows us to clarify the role of different components of the contraction-dominated motility machinery and link 
the proposed mechanism with the previous work on optimization of the crawling stroke irrespective of the underlying 
microscopic processes (e.g. DeSimone and Tatone, 2012; Noselli et al., 2013). It does not, however, address directly 
the main issues of motility initiation and motility arrest that require more elaborate constructions. 

Recall that in crawling cells, the ‘motor part’ containing contracting cytoskeleton (lamellipodium), is a thin active 
layer located close to the leading edge of the cell, see Fig. 1. We assume that all mechanical action originates in 
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lamellipodium and that from the mechanical viewpoint the rest of the cell, including the nucleus, can be interpreted as 
cargo. The main task will be to develop a model of freely moving lamellipodium which we schematize as a segment of 
active gel in viscous contact with a rigid background. The actin network inside the gel is contracted by myosin motors 
which leads to an internal flow opposed by the viscous interaction with the background. The unidirectional motion in 
a layer adjacent to the background that ultimately propels the cell is a result of the asymmetry of contraction. 



cargo 


motor pail 


Figure 1: Conceptual discrete model of the motility mechanism in a crawling keratocyte cell. 


A toy model elucidating this point involves three rigid blocks of size 4 placed in a frictional contact with a rigid 
support, characterized by the viscous drag coefficient The neighboring blocks are connected by active pullers 
(force dipoles) exerting contractile forces. The essential long range interactions representing global volume constraint 
(due to passive elastic structures and osmotic effects, see Section 3) are modeled by a linear spring with stiffness 
k connecting the first and the last block. To regularize the problem we place in parallel with contractile elements 
additional dashpots characterized by the viscosity coefficient rj. In the absence of inertia, we can then write the force 
balance equations in the form 

-h^xi+k^-^^+Xi-ri^^=0 

-ib^x2 -x\ + X2 - ^ ( 1 ) 

-1,^X2-k^-^-X2-n^^=o, 

where xiif), X 2 {t), x^if) are the current positions of the blocks and Lq is the reference length of a linear spring. This 
spring describes the membrane-cortex ‘bag’ around the lamellipodium allowing the inhomogeneous contraction to be 
transformed into a protruding force. We assume that polarization has already taken place and therefore the contractile 
force dipoles;^! > 0 and ;^2 ^ 0 acting between the two pairs of blocks are not the same;^i ^ xi- The polarization 
itself requires additional constructs and will be addressed later. 

System (1) can be rewritten as three decoupled equations for the length of our active segment L{t) = X'iit) - vi(0, 
its geometric center Git) = (^ 3(0 -f vi(0)/2 and the position of a central block X 2 it) representing the internal flow: 

-4^(1 -I- iyi])L =xi +X2 + 2kiLILQ - 1) 

2h^(l+3ll/ll)G=Xi-X2 ( 2 ) 

- 4^(1 + 3lllll)X2 =Xl -X2 

where Iq = yfnJ^ is the hydrodynamic length scale which will ultimately play the role of a regularizing parameter. 
The first equation shows that the length is converging to a steady state value: 

Loo =Lo[l- (xi +T2)/(2k)]. 


Notice that in order to avoid the collapse of the layer due to contraction, it is necessary to ensure that the spring has 
sufficiently large stiffness k > (xi +X2)I2. We also observe that independently of the value of the evolving length L(0, 
the velocity of the geometrical center of our train of blocks V is always the same 


V = G = 


X\ -X2 

24^(1 + 3 / 2 // 2 )- 


( 3 ) 


One can see that the system can move as a whole only if;^i ^ X 2 , which emphasizes the crucial role for motility of 
the polarization and the ensuing inhomogeneity of contraction. 
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We observe that the middle block moves in the direction opposite to the motion of the center of the system with 
a constant velocity X 2 = -2V. Therefore, it takes a finite time ~ LodOV) for the central block to collide with the 




Figure 2 : (a) Schematic representation of the motion of individual particles (blocks) forming the motor part of a crawler in a steady state regime 
(three particle case). Trajectories in space time coordinates of the particles x\ (magenta, OBCEF), X2 (green, ABDEG) and X3 (red, ACDEG); 
dashed lines show the jump parts of the crawling cycle. Continuous flows have to overcome friction while the jumps are assumed to be friction 
free, (b) A closed loop constituting one full stroke in the parameter space (x2 - xi, X3 - X2). The time of one full stroke (A to G) is Tg = LooIV and 
the distance traveled by the crawler per stroke is VTg = Too. 


block at the rear. Beyond this time, system (1) formally collapses and additional assumptions are needed to extend 
the dynamics beyond the collision point. The origin of the problem is our focus on the layer adjacent to the rigid 
background and the neglect of the global fiow of actin. 

To make the model more adequate we have to take into consideration that while the flow of F-actin (polymer¬ 
ized or filamentous actin) is continuous along the contact surface, the cytoskeletal network disintegrates into G-actin 
(unpolymerized monomers) at the trailing edge and reintegrates from the available G-actin at the leading edge. The 
polymerization induced depletion of G-actin at the leading edge is compensated by the diffusive counter-flow of actin 
monomers from the back of the cell to its front. This counter-flow cannot be described directly in the ID setting. 

It can be modeled, however, in an indirect way by mass and momentum preserving periodic boundary conditions 
allowing F-actin to disappear at the rear and reappear in the front. This situation is rather typical for continuum 
mechanics where unresolved spatial scales are often replaced by balance-law-preserving jump/singularity conditions 
as in the case of shock waves, crack tips and boundary layers. 

More specifically, to account for global circulation (turnover) of the cytoskeleton in a one dimensional setting, we 
assume that there is a singular source of actin at the front of the cell that is compensated by the equivalent singular 
mass sink of actin at the rear of the cell. This assumption allows us to close the treadmilling cycle, even though the 
details of the discontinuous part of the cycle, involving both the polymerization reaction and the diffusive transport 
of monomers, are not explicitly resolved in the model. We essentially postulate that there is a pool of G-actin which 
is replenished as fast as it is depleted and that the resulting reverse fiow of actin is synchronous with the direct 
fiow. Under these assumptions the reverse fiow is viewed as passive and and is assumed to be driven exclusively 
by inhomogeneous contraction. In particular, we neglect active propulsion on free boundaries due to growth and 
lessening of the network. 

We describe these processes in our toy model by assuming the possibility of creation and destruction of the blocks. 
Our goal is to account for the fact that actin polymerizes at the leading of the cell (where blocks are assembled) and 
depolymerizes at the trailing edge of the cell (where blocks are disassembled). We offer two interpretations of the 
underlying continuous process in terms of discrete blocks emphasizing different sides of such passive treadmilling. 

In a first interpretation, we assume that as a result of each collision a block at the rear is instantaneously removed 
from the chain and at the same time an identical block is added at the front. In other words, each (equilibrium) de¬ 
polymerization event at the rear is matched by an (equilibrium) polymerization event at the front. Here we implicitly 
refer to the existence of a stationary gradient of chemical potential of actin monomers and of a large pool of monomers 
ready to be added to the network at the front as soon as one of them is released at the rear. The ’instantaneous’ 
reappearance of the disappearing blocks should be understood as a mean to model the overall continuity of the fiow. 
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The structure of the resulting stroke in the t, x plane and in the X2 - xi, X3 - X2 plane is shown in Fig. 2 . One 
can see that each block maintains its identity through the whole cycle and that its trajectory involves a succession 
of continuous segments described by ( 1 ) that are interrupted by instantaneous frictionless jumps from the rear to the 
front. Notice that in this interpretation the blocks can change order and the condition xi < X 2 < X 3 is not always 
satisfied. For instance, when the blocks x\ and X2 collide at point B, the block xi disappears at the back (point B) and 
reappears at the front (point C) ahead of the block V 3 . This jump mimics the frictionless part of the treadmilling cycle. 
Similarly, when when the block V 3 collides with the block X 2 at point D, the latter reappears at the point E ahead 
of the block x\. This interpretation is attractive because it allows one to trace the trajectories of the blocks through 
subsequent treadmilling cycles. It is, however, a bit misleading because in reality the block that disappears at the back 
and the block which instantaneously reappears at the front are definitely not the same even though they are identical. 

According to a second interpretation, illustrated in Fig. 3, the middle block is the only one to undergo cycling 
motion. As a result, the ordering x\ < X 2 < X 3 is always preserved and the distances between the first two blocks 
h = X 2 - xi and the last two blocks 12 = ^ 3 - X 2 can be only positive. We can alternatively say that now the notations 
xi,X 2 , X 3 indicate positions only and can refer to different blocks in different times. In this interpretation, when the 
middle blocks hits the rear one, it is the middle block that gets recycled to the front while the rear one keeps moving 
continuously. 



(a) 


without friction 
-2V 

^ with friction ^ 


(b) 


Figure 3 : Schematic representation of the continuous {a,p) and the jump (/ 5 , cr) part of the crawling stroke. The zero area loop in li,l2 plane 
illustrating the stroke is shown in (a). The loop is not symmetric because continuous flow have to overcome friction while the jumps are assumed 
to be friction free. The ‘tank thread’ analogy in (b) is not fully adequate because the ‘departing’ blocks at point p, that enter the pool of actin 
monomers, and the arriving blocks at point cr, that are simultaneously taken from the same pool, are not the same. 


In coordinates {h^h) the cycle collapses on a single line, which is traveled continuously in one direction and 
discontinuously in the other direction, see Fig. 3(a). Notice that the internal parameters li{t) and l2it) undergo a 
periodic sequence of extensions and contractions which resemble the mechanism propelling the swimming sheet 
(Taylor, 1951) and its crawling analogue (DeSimone and Tatone, 2012). The main difference is that in our case the 
propulsion is achieved because of the asymmetry of friction forces acting in the different phases of the stroke. More 
specifically, we assume that during the continuous phase of the cycle the blocks move with friction (polymerized 
filaments experience effective drag transmitted by focal contacts), while during the discontinuous phase the dissipation 
(associated with reaction and diffusion) can be neglected. The situation is remotely analogous to that of a rotating 
‘tank tread’, see Fig. 3(b), even though in reality the disappearing block and the appearing block are not the same. 
This interpretation is closer to the physical picture where the points of the membrane (cortex) represented by two 
side blocks move with a constant speed ensuring that the cell maintains its length. We reiterate that both discrete 
interpretations are schematic and will be backed later in the paper by an appropriate continuum modeling. 

Since the obtained expression for velocity (3) remains finite in the limit /o/4 ^ 0 it appears that the dashpots play 
a redundant role in this model and can be dropped. To illustrate the role of the dashpots we now consider the case of 
N coupled blocks. Then, the force balance for the central blocks 7 e [2, A - 1] reads 

Xj - Xj-\ Xj - Xj+i 

-Ib^Xj -Xj-I +Xj - n - -J - ri - - -= 0. 

Lb Lb 


1 









This system of equations can be written in the matrix form, 


Tx = b, 


(4) 


where we denoted by x the unknown vector The tri-diagonal matrix 


-( 2 + 1 ) 1 0 0 

‘'0 

1 -( 2 +|) 1 0 

‘‘0 


0 

0 


0 

0 

0 


0 

0 1 -( 2 +|) 1 

0 0 1 -( 2 + 1 ) 

^0 


describes the viscous coupling and frictional interaction with the background while the vector 


-Ti + ^^0 - fxi 
Xi -X2 


b = 


4 

^^0 


Xn-2 -Xn-1 
^/2 

Xn-1 - o-q - -j^xn 


with (To = -k(xN - xi - Lo)ILo carries the information about the active forcing, the mean field type elasticity and 
the boundary layer effects. To find the solution x, we need to invert the matrix T and then solve a system of two 
coupled linear equations xi = (R b)i and % = (R b)A) where R = T“4 The components of the matrix R can be found 
explicitly (Meurant, 1992) 


^ _ cosh i)A) - cosh ((A^ -r 1 - I 7 - /|)T) 

2 sinh(/l) sinh((A^-I-1)/1) 

where A = arccosh(l -r /^/(2 /q)). Knowing the ‘velocity field’, we can now compute the steady state value of the length 

( Zfj/cosh(A(J-N/2))xj) 

Loo — Lq I 1 — ^ I . 

I Zf=i‘cosh(^(y-yV/2))A: j 

From this formula we see again that a finite stiffness is necessary to prevent the collapse of the system under the 
action of contractile stresses: assuming for instance that Xi = T we obtain the low bound for the admissible elasticity 
modulus k >x. 

The steady velocity V = (x^ + x\)12 of the geometrical center of the system can be also computed explicitly 

277 sinh(/lA^/2) 

When N is even, by denoting M = A^/2, we can rewrite this expression in the form 

h iJfji sinh(7A)CtM+i -XM-j) 

V = - - -. 

27] sinh(/lM) 
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from where it is clear that (as in the case of three blocks) the symmetry of the vector;^ with respect to the center must 
be broken for the system to be able to self-propel. 

If we now formally drop the dashpots by assuming that /q = 0 we obtain similar expressions for the velocity and 
for the steady state length as in the three block (A^ = 3) case 


y = 


Xn-\ -Xi 


,L 


OO — 



Xi -^Xn-i \ 
2 k )' 


(5) 


The reason behind this similarity is that, in this limit, the ‘flow’ fully localizes in the two boundary elements, the only 
ones present in the case N = 3. More precisely, the solution of the discrete problem depends singularly on the ratio 
/^//q and becomes progressively more concentrated around the boundary elements as /^//q ^ oo. Such localization 
presents a certain analytical problem if we consider the continuum limit when N ^ oo and 4^0 while Nl^ L, 
where L is the continuum length of the self-propelling segment. Indeed, in this limit the size of boundary layers tends 
to zero and the discrete solution converges to a distribution. The viscosity, introducing a length scale 4, is thus needed 
to preserve the regularity of solutions in the continuum limit. 

Observe also that the limits 4^0 (dropping dashpots) and 4^0 (continuum approximation) do not commute. 
For instance, if we choose in (5) the motor distribution with dllxi = 0 except for one ;^2 = X* > 0 we obtain V = 0 for 
any value of 4, in particular, when 4 ^ 0 we still have V ^ 0. If instead we first perform the continuum limit while 
keeping 4 finite we obtain 


fXcosh[(x-L^/2)/lo]xix)dx 
2 klo sinh [Loo/( 2 /o)] 


( 6 ) 


and 

C “ sinh [(x - Loo/2)/lo]x(x)dx 

V = -- -. (7) 

2t] sinh [Loo/( 24 )] 

If we now take a distribution of motors x(^) = X*^o where Sq is the Dirac mass at v = 0, which can be viewed as a 
continuum analog of the discrete distribution considered above, we obtain that V = ;^*/(2/q^). Then in the limit 4^0 
we obtain that V ^ oo which is in conflict with our previous assertion that V = 0, obtained when the order of limits 
was reversed. Assume now that 4 ~ and hence /^//q ~ 1I(tjN^). One can see that the crossover scaling rj ~ N~^ 
separates the two non-commuting limiting regimes. For /^//^ ^ oo (which is a dimensionless version of 77 N~^) the 

internal flow localizes in the boundary layers whose thickness disappears when 77 ^ 0 ; when we dropped the dashpots 
in the three element model we could not detect this localization because the two boundary links were the only ones 
present in the system. In the other limit /^//q ^ 0 (dimensionless version of 77 » N~^) the viscosity dominates the 
dynamics and the internal flow becomes uniform. 

In the next sections the formulas ( 6 ) and (7) will be obtained directly from the continuum model. We will also 
see more clearly how the introduction of the viscosity-related internal length scale and the associated nonlocality 
regularizes the continuum model which otherwise has only singular solutions. 


3. The continuum model 

We model the lamellipodium as a one dimensional continuum layer in frictional contact with a rigid background, 
see Fig. 4. Assuming that the drag is viscous and neglecting inertia we can write the force balance in the form 

d^o- = ^v, (8) 

where cr(x, t) is the axial stress and v(x, t) is the velocity of the cytoskeleton (actin network). Eq. ( 8 ) is the continuous 
analog of the system (4) in the discrete problem. 

As in the discrete model, we denoted by ^ the coefficient of viscous drag. Such representation of active adhesion 
is usual in the context of cell motility (Rubinstein et al., 2009; Larripa and Mogilner, 2006; Julicher et al., 2007; 
Shao et al., 2010; Doubrovinski and Kruse, 2011; Hawkins et al., 2011). It implies that the time-averaged shear stress 
generated by constantly engaging and disengaging focal adhesions is proportional to the velocity of the retrograde 
flow, see Tawada and Sekimoto (1991) for a microscopic justification. There is evidence (both experimental (Gardel 
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Figure 4 : Schematic representation of a continuum model simulating lamellopodial contraction-driven crawling. 


et al., 2008, 2010; Mogilner, 2009; Bois et al., 2011; Schwarz and Gardel, 2012) and theoretical (DiMilla et ah, 1991; 
Mi et ah, 2007)) that this assumption describes the behavior of focal adhesions accurately only when the retrograde 
flow is sufficiently slow. The behavior of adhesion strength in the broader range of velocities is biphasic and since we 
neglect this effect, we potentially misrepresent sufficiently fast dynamics. Observe though that for both keratocytes 
and PtKl cells the rate of lamellar actin retrograde flow varies from 5 to 30nm.s“^ in usual experimental conditions 
(Schwarz and Gardel, 2012) and in this range a direct proportionality relationship between traction stress and actin 
retrograde flow has been confirmed experimentally (Gardel et al., 2008; Fournier et al., 2010; Barnhart et al., 2011). 
The characteristic velocity of the flow in our problem is 20nm.s“^ which falls well into the aforementioned interval 
where the biphasic behavior can be neglected. 

Following Kruse et al. (2006); Julicher et al. (2007); Bois et al. (2011) and Howard et al. (2011), we assume that 
the cytoskeleton is a viscous gel with active pre-stress. We neglect the bulk elastic stresses that relax over a time scale 
of 1 - 10 s (Rubinstein et al., 2009; Wottawah et al., 2005; Kole et al., 2005; Panorchan et al., 2006; Mofrad, 2009; 
Recho and Truskinovsky, 2013) which is much shorter than characteristic time scale of motility experiments (hours). 
We can then describe the constitutive behavior of the gel in the form 

cr = r}dxV+xc, (9) 

where r] is the bulk viscosity, c(x, t) is the mass concentration of motors and ;^ > 0 is a contractile pre-stress (per 
motor) representing internal activity. The constitutive relation (9) generalizes the parallel bundling of dashpots and 
contractile units in the discrete model. The important new element is that the strength of the contractile elements may 
now vary in both space and time. 

In the discrete model the concentration of motors c was a given as a function of v. To obtain a more self consistent 
description we assume that the function c(v, t) satisfies a convection-diffusion equation (Rubinstein et al., 2009; Bois 
et al., 2011; Barnhart et al., 2011; Wolgemuth et al., 2011; Hawkins et al., 2011) 


d,c + d^{cv) = Dd^^c, (10) 

where D is the diffusion coefficient. Behind (10) is the assumption that myosin motors, actively cross-linking the actin 
network, are advected by the network flow and can also diffuse which accounts for thermal fluctuations. 

To justify this model, consider a simple mixture model with two species representing attached and detached 
motors. The attached motors are advecting with the velocity of actin filaments and can detach. The detached motors 
are freely diffusing, and can also attach. Suppose that the attachment-detachment process can be described by a first 
order kinetic equation. Then the system of equations governing the evolution of the concentrations of attached Ca and 
detached cj motors can be written as: 


dtCa + dxiCaV) = ^onQ “ 


where ^on and are the chemical rates of attachment and detachment and D is the diffusion coefficient of detached 
motors in the cytosol. Now suppose that the attachment-detachment process is chemically equilibrated and hence 
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CalCd = K, where K = is the reaction constant. Then for the attached motors performing contraction we 

obtain 

-^^dtCa + dxiCaV) - —d^xCa = 0 . 

K K 

Our equation (10) is obtained in the limit K ^ oo (fast attachment) and D/K ^ D (fast diffusion). 

Denote by L(t) and l+(t) the rear and front boundaries of our gel segment. To account for cortex/membrane 
elasticity we assume, as in the discrete model, that the boundaries are linked through a linear spring (Barnhart et al., 
2010; Du et al., 2012; Loosley and Tang, 2012; Recho and Truskinovsky, 2013). This assumption affects the values 
of the stress in the moving points L(t) and l+(t): 

(T(U(t),t) = -k(L(t)-Lo)ILo. 


Here L(t) = l+(t) - L(t) is the length of the segment, k is the effective elastic stiffness and Lq is the reference length. 
As we have seen in the discrete model, the presence of an elastic interaction plays a crucial role in preventing the 
collapse of the segment due to contractile activity of motors. 

Our next assumption is that the external boundaries of the self propelling segment are isolated in the sense that 
they move with the internal flow /+ = v(/+). We imply here that the addition and deletion of F-actin particles inserted 
at the front and taken away at the rear does not contribute to propulsion. We also impose a zero exterior flux condition 
for motors dxc(l+(t), t) = 0 ensuring that the average concentration of motors 



c(x, t)dx 


( 11 ) 


is conserved. To complete the setting of the problem we need to impose the initial conditions /+(0) = /+ and c(x, 0) = 
c^(x). 

Our assumption that the bulk stiffness of the cytoskeleton is equal to zero (also known as the infinite compress¬ 
ibility assumption (Julicher et al., 2007; Rubinstein et al., 2009)) allowed us to uncouple the force balance problem 
(which becomes statically determinate) from the mass transport problem. Then by solving the main system of gov¬ 
erning equations (8)-(ll) we can obtain the velocity field and the concentration of motors. To recover the mass 
distribution of the cytoskeleton we need to solve a decoupled mass balance equation with a kinematically prescribed 
velocity field (Recho and Truskinovsky, 2013; Recho et al., 2013). 

Indeed, suppose that by solving the system (8)-(ll) we found the velocity field v(v, t). This means that we also 
know the trajectories of the free boundaries L(t) and l+(t). To find the mass density of actin p(x, t) in the gel, we need 
to solve the mass balance equation 

dtp + dxipv) = 0 ( 12 ) 

with initial condition p(x, 0) = Pq{x). Here we neglected the diffusion of actin which is weak comparing to the 
diffusion of myosin. Now, since both the leading and the trailing edges of the moving lamellipodium coincide with 
the trajectories of particles, the total mass M is conserved 

rhit) 

M = I p(x, t)dx. 

Jut) 

To address the problem of continuous circulation and to close the cycle of the cytoskeleton flow we need to interpret 
the points of density singularities as actin (mass) sources and sinks. In Section 7 we show how the solutions can be 
regularized if we cut out small regularization domains around the sources and sinks and appropriately reconnect the 
incoming and the outgoing flows of matter. 

Dimensionless problem. If we now normalize length by Lq, time by stress by k, concentration by cq and density 

by MILq, we can rewrite the main system of equations in dimensionless form (without changing the notations) 


-Zdxxcr -\-cr = Pc, 
dfC -r n<dx{cdx(r) = dxxC, 


(13) 
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Here we introduced three main dimensionless constants of the problem: 

Z = riK^Ll), 

the ratio of viscous and elastic length scales; 

the ratio of stiffness induced agglomeration over diffusion and finally 

r = cox/k, 


the dimensionless measure of motor contractility. One can discern in (13) the structure of the Keller-Segel system 
from the theory of chemotaxis (e.g. Perthame, 2008). The role of the distributed chemical attractant is played in our 
case by the stress field cr whose gradient is the driving force affecting the ‘colony’ of myosin motors. 

The main mathematical difference between our formulation and the standard chemotaxis problem is that we have 
free boundaries. Using dimensionless variables we can rewrite the boundary conditions in the form 


dxc(l±(t), t) = 0. 


The integral constraint (11) reduces to 

W+ 

J c(x, t)dx = 1. 

In dimensionless variables the mass balance equation (12) takes the form 

dtp + ^dxipdxcr) = 0 , 

and the total mass gets normalized 


rhit) 

t 


p{x, t)dx = 1. 


(14) 

(15) 

(16) 

(17) 


(18) 


Non local reformulation. Since the first of the equations (13) is linear, it can be solved explicitly for cr 


cr(x, t) = 


(L-l)cosh[(G-x)/VZ] P 


cosh[L/(2 VZ)] 


r 


+ — I 'V{x,y)c(y)dy, 


(19) 


where 


^ sinh[(/+-x)/VZ]sinh[(y-L)/VZ] 

'P(x,y) = -—- H{y - x) sinh[(y - x)/ VZ]. 

sinh(L/ yX) 

We introduced the notations: H{x) - the Heaviside function and G{t) = + /+(0]/2 is the position of the geometric 

center of the moving fragment. By eliminating cr from (13)2 we obtain a single non local partial differential equation 
with quadratic non linearity for c(x, t) 

dAx, t) - <K{L - \)dAe{x)c{x, t)] + ^(X, y)c{y, t)c(x, t)dy) = 5,,c(x, t), (20) 

where the auxiliary velocity field 

^ sinh[(x-G)/VZ] 

^ cosh[L/(2 VZ)] 

describes advective flow induced by the elastic coupling between the rear and the front of the active segment. The 
feedback behind contraction-driven motility is contained in the kernel 


(f(x,y) = - 


cosh[(/^ - X)/ VZ] sinh[(y - /_)/ VZ] 
sinh(L/ VZ) 


+ Hiy - x) cosh[()' - x)/ VZ], 


which is due to viscosity-induced bulk interactions in the system and the effect of the boundaries. Notice that this 
kernel has the action/reaction symmetry ^(x, y) = -^(/+ + /_ - x, /+ -r /_ -y) which is a fundamental constraint imposed 
by the balance of momentum (Kruse and Julicher, 2003; Kruse and Julicher, 2000; Torres et al., 2010). 
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Inviscid limit. To distinguish the bulk mechanical interactions from the effects of the boundaries, we use the following 
asymptotic expansion (Ren and Truskinovsky, 2000) 

^b(y - ^) = lim ip{x + G, y + G) 

L—>oo 

= i/ exp(^)ifx-y<0 . (21) 

“ 2 I _exp(^)ifx-y > 0 

In Fig. 5, we compare our viscosity induced interaction kernel with a long range kernel proposed in Kruse and Jiilicher 



Figure 5: Comparison of the bulk part of the viscosity induced interaction kernel ipi, (continuous line) with its mean field analog ips (dashed line) 
proposed in Kruse and Julicher (2003). 


(2003); Kruse and Julicher (2000); Torres et al. (2010), as a model of steric interactions between actin filaments with 
half size Is 


-y) = 


^sgn(y-x), if\x-y\ < Is 
0, if \x -y\ > Is 


( 22 ) 


The length Is plays in Kruse and Julicher (2003); Kruse and Julicher (2000); Torres et al. (2010), the same role as our 
viscous length Iq = represented in (21) by the dimensionless parameter 

Recall that in Section 2 we anticipated a non trivial limit in the continuum theory when the bulk viscosity rj goes to 
zero. Now we see that when Z ^ 0 the kernel becomes singular and the nonlocality in the mechanical part of the 
model disappears. From (13) we also notice that parameter Z enters as a coefficient in front of the highest derivative. 
Therefore, outside the boundary layers of size ~ Vz one can formally assume that cr = Pc which makes the main 
dynamic equation (20) local 

dtc(x, t) + Pn<d^{c{x, t)d^c{x, 0) = dxxc{x, t). (23) 


At small Z the non-bulk contributions to the kernel (p{x, y) will play a role only around the extremities of the moving 
segment and in the limit Z ^ ^ will affect only the boundary conditions. 

By using a new variable w = 1 - we can rewrite Eq. (23) in the form 


dtw(x, t) + djc(wdjcw(x, t)) = 0. 


Here we recognize the porous flow equation which is, however, unusual because the field w(x,t) may be sign- 
indefinite. In particular, in the regimes with c > {%P)~^ one can expect an uphill diffusion similar to that of spinodal 
decomposition. A systematic study of the inviscid case, requiring the knowledge of the boundary conditions in the 
limiting problem, will be done elsewhere. 

Cell velocity. Using the boundary conditions (14) we find from (19) an explicit formula for the (time dependent) 
velocity of the center of our active segment (see also Eq. (7) in Section 2) 


G 


'K'P r‘* 

2Z J,_ 


sinh(^) 

- - — c(x, t)dx. 

sinh(^) 


(24) 
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Similarly we obtain an equation for the evolving length of the segment (see also Eq. ( 6 ) in Section 2) 

nr T nrrp ri+ cosh(^^) 

i = -2—(L-l)tanh(—=)-— c(xj)dx. (25) 

VZ 2^/Z Z Jl cosh(^) 

Notice that in (24) only the odd component of the function c(x, t) (with respect to the moving center G(t)) contributes 
to the integral while in (25) only the even component matters. In particular, if the concentration of motors is an even 
function of v then G = 0 and the segment does not move as a whole. This conclusion is a direct analog of Purcell’s 
theorem (Purcell, 1977) in the case of a crawling body. Notice that for crawling the emphasis is made on spatial 
asymmetry which replaces the emphasis on temporal asymmetry in Purcell’s interpretation of swimming. 

From (24) we infer that the maximal velocity of the self propelling segment is equal to 7C!P/(2Z). If we use 
the data from Julicher et al. (2007); Bois et al. (2011); Howard et al. (2011), we obtain the estimates ;^co - 10^ Pa, 
Lo = 10 yum and 77 = 3 x 10"^ Pa s. For the maximal velocity, this gives xLqCqI(I rj) ^10 yum/min which is rather 
realistic in view of the data presented by Jilkine and Edelstein-Keshet (201 1) and Schreiber et al. (2010). 


Traveling waves. Given our interest in the steady modes of cell motility, which are typical for keratocytes (Barnhart 
et al., 2011), we need to study the traveling wave (TW) solutions of the main system (13). To find such solutions we 
assume that the front and the rear of our segment travel with the same speed i±(t) = V, ensuring the constancy of 
the length L(t) = L, and that both the stress and the myosin concentration depend on x and t through a combination 
u = (x- Vt)/L only. Using this ansatz we find that the equation (13)2 can be solved explicitly 


^ ^ exp[5’(M) - VLu] 

c(m) = —i- 

L Jq exp[s(u) - VLu]du 


(26) 


Here for convenience we introduced a new stress variable s(u) = [cr(u) + (L - 1)] which represents the 
geneous contribution to internal stress field due to active pre-stress. The system (13) reduces to the single 
equation 


~s"{u) + s{u) ■ 


TC{L - 1) = TCP 


exp[5’(i/) - LVu\ 

L Jq exp[5'(t/) - VLu]du 


inhomo¬ 

nonlocal 

(27) 


supplemented by the boundary conditions 


^(0) = ^(1) = 0 and s'(0) = s'(l) = LV. 


(28) 


The two ‘additional’ boundary conditions in (28) allow one to determine parameters V and L along with the function 
s(u). After the problem (27, 28) is solved, the motor concentration profile can be found explicitly by using Eq. (26). 


4. Static regimes 

Initiation of motility is associated with a symmetry breaking instability of a static (non-motile) configuration. To 
identify non-motile configurations we need to find solutions of (27) with V = 0. These solutions may still characterize 
the states with nontrivial active internal rearrangements of both actin and myosin. Static solutions with periodic 
boundary conditions were studied in Bois et al. (2011) and here we complement and extend their analysis. 

If y = 0, Eq. (27) simplifies considerably 

-^s" + s-TC(L-l) = TCP-- (29) 

^ lJ^ exp(s(u))du 

The nonlocal equation (29) was studied extensively in many domains of science from chemotaxis (Senba and Suzuki, 
2000) to turbulence (Caglioti et al., 1992) and gauge theory (Struwe and Tarantello, 1998). In our case, this equation 
where parameter L remains unknown, has to be solved with three boundary conditions 5 ''( 0 ) = 5 '( 0 ) = ^'(l) = 0 because 
the forth boundary condition s'(l) = Ois satisfied automatically. 
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Figure 6: Three families of trivial static solutions L+, L_ and Lq parameterized by P. Solid lines show stable branches while dotted lines correspond 
to unstable branches. Arrows depict the basin of attraction of each branch (see Section 6) 


We begin with the study of the regular solutions of (29). Instead of and !P, it will be convenient to use another 
set of parameters A := 7C(L - 1) < 0 and B \= KVI{L exp[s(ii)]dii) > 0. In terms of parameters (A, B) the problem 
(29) reads 

Z 

-^s" -\-s-A = Bqxp(s) with s'(0) = ^(0) = ^(1) = 0. (30) 

A trivial homogeneous solution of this problem s(u) = 0 exists when A B = 0 which is equivalent in the (!P, 7C) 
parametrization to L = L+ with, 

L± = (1 ± Vl - 4P)/2. (31) 

The sub-branches with longer and shorter lengths L+(!P) and L-{P), respectively, that meet at point a where L-(^) = 
L+(P) are illustrated in Fig. 6. 

To obtain nontrivial static solutions we multiply (30) by 5’', integrate and use the boundary conditions to obtain 
the ‘energy integral ’= W(s), where 

W{s) = ~ “ !])• 

The general solution of this equation can be expressed as a quadrature, u = ± W~^^^(r)dr where we recall u 

designates the normalized space variable. A detailed analysis of this equation is given in Appendix A, where different 
families of static solutions are identified as 5 ^ and (S ^)' where the index + specifies the L+ trivial branch from which 
a particular solution bifurcates: the associated lengths L+ are defined in (31). The integer valued index m corresponds 
to the number of spikes in the configuration s(u). The prime differentiates between two subfamilies belonging to 
the same bifurcated branch with primed subfamily having a length L larger than in the trivial configuration and non- 
primed subfamily having the length L smaller than in the trivial configuration. Fig. 22 illustrates the families S | and 
S 2 - For each family we plot the length of the fragment L as a function of one of the controlling parameters, see Fig. 7. 

In addition to regular solutions described above Eq. (29) has measure-valued solutions corresponding to collapsed 
cells with length Lq = 0- First of all, as we see in Fig. 6, L-(!P) ^ 0 when P ^ 0 (point a') and the limiting 
distribution of motors is concentrated on an infinitely small domain. To characterize the asymptotic structure of the 
singular solutions we suppose that L « 1 and that the maximum of 5' is of order L. Then, by ignoring high order 
terms, we deduce from (29) a simplified boundary value problem 

-s" 7C!PL/(Z X' [1 + S{u)\du) with s'(0) = s(0) = s(l) = 0. (32) 

Then s{u) ^ ^PLu{\ - u)I(IZ)- and the remaining boundary condition ^-^O) = 0 is automatically satisfied in the limit 
L ^ 0. We can then conclude that the singular solutions are of the form 

5’(v) = limL/(v/L) 
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Figure 7: Bifurcation diagram for the L+ branch at fixed V = 0.245 and = 1. See also Fig. 9(a). 


where f{u) = ^Pu(l - u)l(2Z)- Singular solutions of this type can be implicated in the description of cell splitting in 
a cortical geometry (Turlier et ah, 2014); they are also known in other fields where stationary states are described by 
equation (29) (Caglioti et al., 1992; Chen and Lin, 2001; Ohtsuka, 2002; Gladiali et ah, 2012) . The presence of such 
solutions is a sign that in a properly augmented theory, accounting for the vanishing length, one can expect localization 
with active contraction balanced by a regularization mechanism, which may be, for instance, active treadmilling 
(Recho and Truskinovsky, 2013). Our numerical solutions of a non-steady problem, which are naturally regularized 
because of the finite mesh size (see Section 6 ), show that the almost singular solutions of the type described above 
are indeed attractors for initial data with L < L_ when P < 1/4. Moreover, numerical experiments suggest that they 
are the only attractors for P > 1/4. This means that even in the presence of a cortex-type spring, an active segment 
fragment necessarily collapses after the contractility parameter reaches the threshold Pmax = 1 /4. 


5. Re-entrant bifurcations 

We first show that motile branches with V ^ 0 can bifurcate only from trivial static solutions with s(u) = 0, V = 0 
and L = L+. For V 9 ^ 0 equation, multiplying (27) by 5 *' - VL, we find the relation: 

1 - exp(-Ly) = LV r Qxp[s(u) - VLu]du. (33) 

Jo 

Then in the limit V ^ 0 we obtain that Qxp(s(u))du = 1. Since static solutions s(u) must be necessarily sign 
definite, see Appendix A, Eq. 33 implies that the corresponding static solution can only be trivial s(u) = 0. As we 
have seen in Fig. 6 , there are two non-singular families of trivial solutions: one with longer (L+ family) and the other 
with shorter (L_ family) lengths. 

Bifurcation points. To find the bifurcation points along the trivial branch ( 5 - = 0, V = 0, L = L+{P)), we introduce 
infinitesimal perturbations dsiu), dV, SL and linearize (27) together with boundary conditions (28). We obtain the 
boundary value problem 


Ss" - o/Ss = 


Zo?-L^l^2L-l L3(L-1) 


P{L-\) 


Z — —oj^SL -f 


( 2 w- i)dy , 


Here we introduced the notation 


Ss(0) = Ss(l) = 0, Ss'(0) = S s'(1) = LSV. 


0? = (L^ - TCPQIZ. 


(34) 

(35) 

(36) 
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Figure 8: Locus of the bifurcation points in the (TC, V) plane fox X = 1. Insert shows a zoom on the D\ branch around the turning point at ^ = 1/4 
where L_ and L_ branches meet. The detailed bifurcation diagrams for P = 0.245 and = 70 and 100 are shown in Fig. 9(a) and Fig. 10 from 
where the meaning of labels y, y becomes clear. The bifurcation points related to the cut = 2600 (red dashed line) in the {V, L) space are 
shown in Fig. 6. 

Since Cc) = 0 at the trivial branch Ss = SV = SL = 0, we can assume that oj ^ 0. The general solution of the problem 
(34), (35) can be written explicitly 


Ss(u) = Cl smh(-(jiju) + C 2 cosh(-oju) - 


Zo? - D- 
<xP-L\L-\) 


z- 


.2L- 


1 2 .. , l\l- 
—(jj oL + — 


1 ) 


( 2 m - l)SV 


From boundary conditions, non trivial solutions branch from the trivial ones if the matrix 


1 

0 

{2 L-\){D-o?Z) 

(L-1)L3 


(*-■)! 

COSh(6c)) 

sinh(6t;) 

{2L-l){L^-o?Z) 

(7-1)73 

\L\ 


0 

OJ 

0 


_ !l_ 

uP-X 

OJ sinh(6t;) 

OJ cosh(u;) 

0 


73 

/ 


has a zero determinant. This gives a transcendental equation for oj 

2L[cosh(oj) - 1] - ^Poj sinh(6t;) = 0. 


(37) 


(38) 


(39) 


The detailed analysis of this equation is presented in Appendix B. The full locus of bifurcation points in the 
(7C,^) plane is shown in Fig. 8. The lines of bifurcation points + and - originating on the trivial sub-branches 
L+ and L_ smoothly connect at ^ = 1/4, see Fig. 6. When parameter P is held constant while is changing 
each family Df and Si in Fig. 8 is represented by two points. For solutions bifurcating from the trivial branch L+, 
we have 7C+ = (L+ - Xo?)KPL+), which gives points • • • and for the branch L_, we have = 

{t}_ - Xo^^)l(!PL-) which gives points ^^2 ’ • • • Notice that the total number of bifurcation points increases 

to infinity as ^ 00. Now consider the case when = const and P is varied. A line = const in the (7C, P) 
plane cuts again each line of the bifurcation points Di and S i in two points which we denote Z)*, 5 J,... (solutions with 
longer lengths) and 5^*,... (solutions with shorter lengths), see Fig. 6 and Fig. 8. In most cases, one of these two 
points is a bifurcation originating from the L_ trivial solution while the other is from the L+ trivial solution. However, 
as we show in the inset in Fig. 8, the two points may also bifurcate from the same branch L+. As we show later in the 
paper such bifurcations are of particular interest because they describe both motility initiation and motility arrest. 
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Figure 9: (a) Bifurcation diagram with as a parameter showing nontrivial solutions branching from families of homogeneous static solutions 
L+ and L_. The value P = 0.245 and = 1 are fixed. Solid lines show stable motile branches while all the dotted lines correspond to unstable 
solutions. The internal configurations corresponding to branches indicated by numbers (1,1', 2,2', etc.) are shown in Fig. 9(b). The projection 
of the bifurcation diagram on the (TC, L) plane is also shown below, (b) Internal profiles associated with successive bifurcated solutions shown in 
Fig. 9(a) for P = 0.245 and ^ = 1. Our notation (1,3) correspond to asymmetric motile branches while (2,4) describe symmetric static branches. 


Structure of bifurcations. After the bifurcation points are known one can use the Lyapunov-Schmidt reduction tech¬ 
nique to identify the nature of the corresponding bifurcations (Nirenberg, 1974; Koiter, 1976; Amazigo et al., 1970). 
The analysis presented in Appendix C shows that the bifurcations from the trivial to the nontrivial static branch are 
always transcritical. The bifurcations to motile branches can be either subcritical or supercritical. In particular, at 
a given the bifurcation from a static homogeneous solution with longer length is always supercritical while the 
bifurcation from a static homogeneous solution with smaller length can be either subcritical or supercritical depending 
on the value of see Appendix C. 

Bifurcated branches. To illustrate different types of bifurcations we constructed the nonlinear continuation of the 
bifurcated branches by solving the boundary value problem (27)-(28) numerically for successive values of parameters 
and T (tracking algorithm, see Doedel et al. (2007)). In Fig. 9(a), we show the continuation in for both static and 
motile configurations at fixed the corresponding profiles of motor concentration, stress and velocity are shown in 
Fig. 9(b). One can see that each pitchfork (for motile branches) and each transcritical (for static branches) bifurcation 
points gives rise to two nontrivial solutions. For instance, along the static branch L+, the bifurcation point Z)| is 
associated with two motile supercritical branches whereas the point is associated with two transcritical static 
branches. Each pair of motile solutions is symmetric with two opposite polarizations corresponding to two different 
signs of the velocity. Along the first motile branch originating at Z)|, the myosin motors concentrate at the trailing 
edge. For the second motile branch originating at there is an additional peak in the concentration profile, see 
Fig. 9(b). In contrast, the static bifurcation point 51 gives rise to two symmetric configurations with different lengths 
and with myosin motors concentrated either in the middle of the cell or near the boundaries, see Fig. 9(b). As one 
would expect, the higher order static and motile bifurcation points produce solutions with more complex internal 
patterns. For the branches bifurcating from the trivial configurations belonging to L_ family, the picture is similar, see 
Fig. 9(a). 

In Fig. 10, we show in more detail the nontrivial solutions originating from the motile bifurcation points Di at 
two values of parameter which correspond to two sections ajS and a0 shown in Fig. 8 (insert). Notice that a single 
solution connects the bifurcation points D\ (suprecritical) and Z)** (sub- or supercritical) which may belong either to 
one family L+ (aJS where is the same as and is the same is T)[) or to two different families L+ and L_ (a0 
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Figure 10: Bifurcation diagrams along parameter P showing motile branches connecting points and £>**. Corresponding bifurcation points are 
shown in insert in Fig. 8. Solid lines show stable motile branches while all the dotted lines correspond to unstable solutions. The projection of the 
bifurcation diagram on the (P,L) plane is also shown. Parameter is fixed in each graph to = 70 and = 100. Internal profiles on the two 
symmetric motile branches are also shown for = 100. Parameter = 1. 

where D\ is the same as and DJ* is the same as Z)“). In the former case, the nontrivial motile branch has a turning 
point at a finite value of ^ < 1/4 giving rise to a re-entrant behavior. Similar behavior was also observed in some 
other nonlocal models (e.g. Kruse and Julicher, 2003; Tjhung et al., 2012; Giomi and DeSimone, 2014). 

As illustrated in Fig. 10 and shown more clearly in a phase diagram in Fig. 11(a), in the re-entrant regime (suf¬ 
ficiently low 7C), the increase of the average concentration of myosin (increase of P at fixed 7C) first polarizes the 
cell and initiates motility, but then, if the contractility is increased further, the cell may becomes symmetric again by 
re-stabilizing in another static homogeneous configuration (see Fig. 10, = 70). We reiterate that re-symmetrization 

and arrest prior to division (known also as ‘mitotic cell rounding’) is a common feature of almost all animal cells 
(Stewart et al., 2011; Lancaster et al., 2013; Lancaster and Baum, 2014). In this respect, it is interesting that if con¬ 
tractility {P) is increased further, the cell collapses to a point because our effective ‘size preserving spring’ cannot 
support the contraction any more. Following Turlier et al. (2014), we can associate such collapse with cell division. 
We can then argue that our deliberately minimalistic model succeeds in reproducing a rather general pattern of cell 
behavior by showing that symmetrization (stabilization) in space immediately precedes the division. 

While the physical meaning of the non-dimensional parameter !P in this discussion is rather clear (contractility 
measure), the significance of varying at fixed Z. is less obvious because both of these parameters depend on 
frictional strength of the background. Adhesivity of the cell to the substrate is known to be a crucial parameter for 
motility initiation and arrest for various cell types (Banerjee and Marchetti, 2011; Lober et al., 2014). To explicitly 
expose the role of friction, it is instructive to interpret parameter 1 /7C as a measure of adhesivity while maintaining at 
a constant level the parameter which does not have any relation to surface friction and is just a ratio of diffusive 
and viscous time scales. 

The resulting phase diagram in the (^, 1/9C) plane at fixed is shown in Fig. 11(b). In this diagram a 

horizontal path extending from left to right means fixed adhesivity and increasing contractility. One can see that at 
high adhesivity motility regimes cease to exist and static solutions collapse as contractility increases. If the adhesivity 
is below a certain threshold, the contractility increase first causes polarization of a static configuration and motility 
initiation; further increase of contractility causes re-symmetrization, arrest and eventually collapse. An interesting 
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Figure 11: (a) Phase diagram of the system (13) in the parameter plane (TC, P) at fixed Z = I- (h) Phase diagram of the same system (13) in the 
parameter plane (P, 1 /TC) at fixed Z/'K = 0.015. The solid (red) line indicates the motile bifurcation point similar to Fig. 8), while the black 
dashed lines indicate the collapse threshold (Pmax = 1 /4). 


regime corresponds to the very tip of the motile domain shown in Fig. 1 1(b). Near this ‘critical’ point the motility can 
be sustained in a narrow ‘homeostatic window’ of parameters and can be easily arrested by either increase or decrease 
of contractility. 

Very recently new experimental results elucidating motility initiation in fish keratocytes have appeared (Barnhart 
et al., 2015). According to these experiments, at a fixed contractility level (fixed P in our model), the increase of 
surface adhesivity (increase of 1 /7C in our model) promotes static configurations while lowering adhesivity initiates 
motility. As it follows from Fig. 11(b), these observations are in agreement with our predictions. Our model also 
explains another observation made in Barnhart et al. (2015) that at a fixed adhesivity, a blebbstatin (a contractil¬ 
ity inhibitor) treatment promotes arrest of the cells while a calyculin A treatment (a contractility stimulator) initiates 
motility. The question whether a more substantial increase of contractility in experiment can lead to re-symmetrization 
and arrest remains open. It is promising in this respect that some cells are known to undergo static to motile transfor¬ 
mation in response to a decrease in the level of contractility Liu et al. (2010); Hur et al. (2011). The minimal model 
presented in Barnhart et al. (2015) is exactly a 2D version of the one formulated in Recho et al. (2013) and further 
developed in the present paper. While active protrusion and non linear regulation of adhesion were also accounted for 
in Barnhart et al. (2015) to get a realistic cell shape, it is rather remarkable that the fundamental pattern of motility 
initiation (including its dependence on contractility and adhesivity) can be already captured within our much more 
transparent setting, see Fig. 11(b) and Section 8. 

Nonlinear active stress. The fact that the bifurcation leading to polarization and motility initiation is always a super¬ 
critical pitchfork indicates that this model does not allow for metastability resulting in the coexistence of motile and 
non-motile configurations that was observed in other models (e.g. Ziebert and Aranson, 2013; Tjhung et al., 2012; 
Giomi and DeSimone, 2014). To obtain such a coexistence in the present setting, we need to modify our model only 
slightly. The main idea is to consider a more realistic nonlinear dependence of the active stress on motor concentration 
which is linear for small values of c but then saturates after around a threshold c*. More specifically, we rewrite the 
main system of equations in the form 


-Zdxxcr + cr = PO(rc)/r, 
d,c + 'Kdxicdxcr) = dxxC, 

where, following Bois et al. (2011), we choose a particular form of nonlinearity 0(x) = v/(l -l- x) and where we 
introduce the new non dimensional parameter r = cqIc"". 

For simplicity we analyze below only the ‘rigid’ limit when k ^ oo, L ^ Lq while the stress on the boundaries 
-k(LILo - 1) remains finite. Notice that in this limit, which formally means that P ^ 0 and ^ oo, we have to 
re-scale the stress by cqx instead of k and as a consequence in the rest of the section we denote, with some abuse of 
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Figure 12: Bifurcation diagrams in the nononlinear model with fixed length (infinite stiffness) (41) showing the possibility of a switch from 
supercritical to subcritical bifurcation. Parameters: ^ = 1. 


notations, cr := o-jV. The new dimensionless parameter replacing and ^ is T = cqxK^D) = that is assumed to 
be finite (see also Bois et ah, 2011; Howard et ah, 2011; Hawkins et ah, 2009, 2011). In dimensionless variables the 
residual stress can be written as (Tq = - lim^^o lim/^^i(L - 1)/^. Then the boundary conditions read 


U-L 

t) 

L 


= 0 
= CTo 
= 0 

= AdxCr(U(t),t). 


For TW solutions we can write the analogue of (27) 


-Z/ 


+ iS" + — 



r 


( 


f 

Jo 


exp(5’ - Vu) 
exp(5' - Vu)du ^ 


(41) 


where 5- = A{cr - ctq) and sq = Actq. The boundary conditions take the form 5'(0) = s(l) = 0 and ^s-^O) = ^-^l) = V. 
The difference with our static solutions, described in Section 4, is that now we have to find the stress at the boundary 
So instead of the length L. 

The analysis of the motility initiation bifurcation in this case is presented in Appendix D. The results are illustrated 
in Fig. 12. As we see, when the nondimensional parameter r is small, which means that we are in the linear regime, 
the bifurcation from static to motile regime is a supercritical pitchfork. However, at larger values of r the nature 
of the bifurcation changes from supercritical to subcritical. This creates a domain of parameters where static and 
motile regime can coexist and where the system may exhibit metastability and hysteresis. Another important effect 
is that in this range of parameters the motility initiation/arrest is a discontinuous transition which may explain why 
experimenters were unable to observe particularly small velocities of self propulsion in keratocytes (Barnhart et al., 
2011). An alternative explanation of this experimental fact based on the idea of optimality and compatible with the 
supercritical nature of the motility initiation bifurcation, was proposed by Recho et al. (2014) 


6. Stability of post-bifurcational regimes 


Stability of various branches of the TW solutions identified in the previous sections was studied numerically. 
Since we have to deal with a moving segment, it is convenient to map system (13) onto the fixed domain [0,1] 
which makes the coefficients of the governing equations time dependent. To this end, we introduce the new space 
variable u(x,t) = [x - L(t)]/L(t) e [0,1] and denote the new unknown functions &{u,t) = (t[/_ + L(t)u,t] and 
c(u, t) = L{t)c{l- + L{t)u, t]. Then the original problem (13), (15)-(16) takes the form 


-j 2 ^uuO- + cr= -c 


and d,c + ^5„(vc) = 


(42) 
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Figure 13: Cell length L(r), velocity G{t) and profiles clL{u,t), slL(u,t) and v(u,t) - G(t) for the test with initial data shown at r = 0 with 
L(0) = 0.4. Parameters P = 0.245, TC = 150 and = 1 as in Fig. 9(a). The layer collapses due to the contractile stress. 


Here we defined the relative velocity v := - G - (u- 1/2)L, where 

G = (7C/L) [du&(h t) + 01/2, 

L = i^lL) [du&(h t) - 01. ^ ^ 

The remaining boundary conditions can be written as 

d-(M, t) = -{L - 1) and duc(u, t) = 0 at u = {0,1}. (44) 

while the initial data take the form c(u, 0) = c^(u), G(0) = and L(0) = LP. 

We integrated the dynamical system (42)-(44) with initial data chosen close to one of the known steady states. 
The numerical scheme was based on the finite volume method (LeVeque, 2002). We used two dual regularly-spaced 
grids on the interval [0,1]: Z and Zj. Given the initial condition c we solved (42) i on Z and computed the effective 
drift term v on Zj. We then applied the upwind finite volume scheme to (42)2 and updated the concentration profile 
c on Z which provided us with the new initial data for the next time step. The time interval for each time step was 
adapted to ensure that the Courant-Friedrichs-Lewy condition is uniformly satisfied on Zj. 

Our numerical experiments suggest that the trivial branch L_ is unstable together with all nontrivial non-singular 
static solutions. The singular static solutions from the Lq family appear to be locally stable. To illustrate the attractive 
nature of the singular static solutions we choose in Fig. 13 the initial configuration with a length smaller than L_ with 
an internal initial profile biased to the front associated to a motile solution. We observe that the length collapses to 
zero in finite time and cell velocity goes to zero. In accordance with the computations made in Section 4, the stress 
profile converges to s(u)IL ~ ^Pu(u - l)/2, velocity to v(u) ~ - 1/2) and concentration to c(u) = 1. 

Next, we observed numerically that the dynamic solutions are all unstable except for the branches bifurcating from 
the points on the trivial branch L +. The trivial branch L+ branch is locally stable until the first (motile) bifurcation 
Both symmetric subbranches of (subfamilies 1 and l' in Fig. 9(a) and Fig. 9(b)) are stable. To illustrate the 
instability of a nontrivial static solution, we show in Fig. 15 the escape of the phase trajectory from the neighborhood 
of the trivial static solution L_. Since in this numerical test the value of was chosen to be smaller than the critical 
value, corresponding to the bifurcation of the first motile branch D^, the system originally placed near L_ becomes 
unstable and then re-equilibrates on another trivial static branch L+ without moving its geometrical center. 

In Fig. 14 we illustrate motility initiation in two initially almost identical and nearly homogeneous static configu¬ 
rations which differ by a localized concentration peak introduced either at the rear or at the front of the cell. We see 
that with time these two initial profiles converge to the different stable motile solutions Di and D[. The initial inho¬ 
mogeneity is remembered and selects the subfamily of the Di solutions with the same bias. As we see, independently 
of the direction of motion the cell recovers its length after a short transient period. 
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Figure 14: Cell length L{t), velocity Git) and profiles c, cr and v(m, t) - G(t) for the test with initial data shown at r = 0 with L(0) = 0.5. Parameters 
P = 0.245, TC = 150 and = 1 as in Fig. 9(a). The layer polarizes to one the motile attractor (depending of the initial bias). 



Figure 15: Cell length L(t), velocity G(t) and profiles c, cr and v(u, t) - G(t) for the test with initial data shown at r = 0 with L(0) = L_ and the 
homogeneous concentration c^iu) = 1/L_. Parameters P = 0.245, 9C = 50 and = 1 as in Fig. 9(a). The layer restabilises from the homogeneous 
branch L_ to L +. 
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Figure 16: Cell length Lit) and velocity G(0 for the test with V = 0.245, TC = 400 and = 1 starting from homogeneous initial state with different 
initial lengths L(0) = 0.6 (left) and L(0) = 0.5 (right). The labels refer to Figs. 9(a), 9(b). The two non trivial static branches bifurcating from 
denoted 2 and 2' on Figs. 9(a) and Figs. 9(b)) have very different kinetic properties. 


As in Bois et al. (2011); Howard et al. (2011); Kruse and Jiilicher (2003), who considered the problem with fixed 
boundaries, we find that some unstable multi-peaked static and dynamic solutions are long living. This behavior is 
reminiscent of the spinodal decomposition in a ID Cahn-Hilliard model where the coarsening process gets critically 
slowed down near multiple saddle points (Carr and Pego, 1989). To illustrate the long transients near the unstable 
solutions we study in Fig. 16 evolution of two initially homogeneous concentration profiles with different initial length. 
We observe that the phase trajectory first approaches the unstable branch from subfamily 2 from Fig. 9(a) and Fig. 9(b) 
before being finally attracted by the stable configuration from the subfamily 1'. Interestingly, the symmetric subfamily 
2' can be also initially approached if we choose slightly different initial data, however, this regime is abandoned much 
faster than the solution from the subfamily 2, see Fig. 16(b). Based on our simulation, we conjecture that the lifespan 
of an unstable branch is linked to the distribution of motors and the states with higher localization of motors on the 
periphery of the cell survive longer than the states where motors are spread near the center of the cell. 

To summarize, we found considerable numerical evidence that in a problem Wiih free boundaries only trivial static 
solutions can be stable and only solutions with monotone profiles can describe configurations of steadily moving cells. 
To confirm these results a more systematic mathematical analysis of stability of the obtained TW solutions is needed. 
Cells with constrained or loaded boundaries may show different stability patterns as it is evidenced by the study of 
a related problem with a periodic boundary conditions (Bois et al., 2011; Howard et al., 2011; Kruse and Jiilicher, 
2003). 


7. Mass transport of actin 


As we have already mentioned, the infinite compressibility assumption allowed us to decouple the force balance 
equation from the mass balance equation. Once the velocity field v(x, t) = t) is known, the latter can be 

solved a posteriori by the method of characteristics. 

Denote the trajectories of the mass particles by x = 0(f, t), where /_(0) < f < /+(0) is the Lagrangian coordinate 
ait = 0 and L(t) < 0(f, t) < l+(t). The characteristic curves can be found from the equations 


#(^. s) 

ds 


v(0(f, s), s). 


(45) 


Along these curves we must have 


Jp(0(^, S), S) 

ds 


-p(0(f, s), s)dM(f>(^, sf s), 


Integration of this equation gives gives an explicit formula for the mass density 


p(0(f, 0, t) = po(f) exp{- r s), s)ds}. 

Jo 
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As we are going to see below, this solution is applicable only outside the singular points describing the sinks and the 
sources. 

Consider a TW solution of (13) which satisfies the boundary conditions l-{t) = Vt and l+(t) = L-\-Vt. Introducing 
the normalized co-moving variable ^ = (cf) - Vt) IL and the normalized Lagrangian variable f = ^/L(0), both in the 
interval [0,1], we obtain that v = v(0) and Eq. (45) reduces to 

0 ^ v(0(^, 0 ) - y 

dt L 

For TW solutions the general formula (46) describing the mass distribution simplifies 

Lp(0(f, 0,0{v(0(f)) - V} = L(0)po(f){v(f) - V}. 

According to (47) the points of the body where v = V are singular because the relative flow there is stagnated. If 
at such point the slope of the function v( 0 ) is negative we obtain a sink of particle trajectories 0 = 7 + (i.e. an attractor 
for particles as r ^ 00) whereas if the slope of the function v( 0 ) is positive, the singular point 0 = 7 - corresponds to a 
source of particle trajectories (an attractor as ^ ^ -00). An important feature of the fiows described by (47) is that it 
takes an infinite time for a mass particle to reach a sink or to leave a source because (v( 0 ) - V)~^ is not integrable in 
the neighborhood of 7 _ and 7 +. 

jy- |v( 0 )-y| 

This implies that mass density infinitely localizes in the singular points (sources and sinks) because Lp( 0 )|v( 0 ) -V\ = 

=0. Then all mass points (corresponding to different values of f) come from the sources where the characteristic 
curves accumulate at large negative times and disappear in the sinks where the characteristic curves accumulate at 
large positive time. 

For the trivial static solutions characterized by the lengths L+, there is no flow (v - V = 0) and the mass density 
does not depend on either space or time. The density profiles for nontrivial static and motile solutions can be illustrated 
near the bifurcation points where the velocity profiles are known explicitly. 

For instance, in the case of the nontrivial static branches introduced in Section 5, we obtain 

dStt, t) ^ ^ 

^^=^sin(u;,0(f,O), (49) 

dt 

where cOc = -Imn. For determinacy, we choose the value of the amplitude g in such a way that the maximum of 
our dimensionless velocity field is equal to one. The approximate value of ^ can be computed in the vicinity of the 
bifurcation point from the amplitude equations presented in Appendix C. In Fig. 17(a) we show sample solutions of 
(49) corresponding to homogeneous initial conditions 0 (^, 0) = ^ for positive and negative values of ^ corresponding 
to the two possible branching directions. The corresponding density profiles are illustrated in Fig. 18 where the passive 
treadmilling cycles are shown by arrows. 

Similarly, for the motile branches we need to solve the characteristic equation 

^ K t) - 1/2)) - 2 sin(a;,/2)l - l|, (50) 

dt ( (jjI cos(6t;c/2) L J J 

where ojc is a solution of the equation (B.2). Both equations can be solved analytically by separation of variables. In 
Fig. 17(b), we show the sample solutions of (50) corresponding to homogeneous initial conditions 0(^,0) = I again 
for the positive and negative values of 

We reiterate that this model is singular because in a one dimensional setting we are obliged to over-schematize the 
treadmilling of actin. 

To recover the circulation aspect of the flow in a one-dimensional setting, we need to regularize the problem 
near the singular points and make the mass flux finite. For instance, we can cut out small regularizing domains of 
size 6 around sinks and sources. In this way we obtain an effective ‘polymerization zone’ around each source r_ = 
{0 e [0,1]/|0 - y_| < e| and an effective ‘depolymerization zone’ around each sink r+ = {^ € [0,1]/|^ - y+| < e|. 
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Figure 17: (a) Trajectories of particles from sources to sinks for the first two static bifurcation points for initially homogeneously distributed set of 
particles, (b) Trajectories of particles from sources to sinks for the first two motile bifurcation points for initially homogeneously distributed set of 
particles. Labels 1, l\ 3,3' and labels 2,2', 4, 4 ' are related to Fig. 9(a) and Fig. 9(b). 



Figure 18: Density profiles for the first two motile and static branches for ^ > 0, the profiles for ^ < 0 are the same; only the treadmilling cycles 
(indicated by black circles) are going in the opposite direction. Labels are related to Fig. 9(a) and Fig. 9(b). Parameter is 6 = 0.01. 
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Figure 19: (a) Comparison of numerics with the experiments performed by Verkhovsky and co-authors in Verkhovsky et al. (1999). Parameter 
values: Z. - 0.0125, ^ = 0.1, TC = 20. Integration is started from an initial cell length of L; = 1.12 with an homogeneous distribution of motors. In 
insets we show some snapshots of the distribution of motors obtained by numerical integration. The homogeneous configuration switches fast to a 
symmetric distribution where motors are relocating at both sides forming a two opposed lamellipods system. This state is long living but unstable 
and application of a transient loading leads to a break of symmetry and the subsequent localization of motors to the trailing edge forming a one 
lamellipod system, (b) Locus of the first motile bifurcation point associated to the homogeneous L+ branch for Z - 0.0125. Red dot shows the 
experimental data for keratocyte = 0.1 and TC = 20 which belongs to the motile regime. Such regime would be spontaneously reached under 
infinitesimal perturbations from a symmetric state but the long living nature of regime 2 (see Fig. 16) makes it necessary to impose a transient 
asymetric perturbation to observe motility in experiments. 


We assume that in the domain r_ the network is constantly assembled from the abundant monomers while in the 
domain r+ it is constantly disassembled so that the pool of monomers is replenished. The ensuing closure of the 
treadmilling cycle is instantaneous (jump process) allowing the monomers to avoid the frictional contact with the 
environment. More precisely, we assume that the jump part of the treadmilling cycle is a passive equilibrium process 
driven exclusively by myosin contraction. The turnover time 

Jdr_ |v(0) - y| 

is now finite and the corresponding density profiles are illustrated in Fig. 18. Notice that the fiow between the neigh¬ 
boring source and sink can be interpreted as a treadmilling cluster. Thus, for the static branch, we have 2m such 
clusters and for the m* motile branch we have 2m - 1 clusters. We reiterate that according to the numerical stability 
analysis conducted in Section 6, the only stable motile solutions as the ones with a single treadmilling cycle extending 
from the leading to the trailing edge. 


8. Experimental verification of the model 

We can now compare the predictions of the model with experiments describing motility initiation in keratocytes. 
For instance, in the experiment of Verkhovsky et al. (1999), a mechanical force was transiently applied via a mi¬ 
cropipette on one side of a keratocyte fragment. Since the data presented in Fig. 5 of Verkhovsky et al. (1999) (and 
reproduced with permission in our Fig. 19) are of one dimensional nature we can directly apply our model after 
adjusting it to account for mechanical loading. 

In order to make quantitative predictions we need to specify the values of parameters relevant for fish keratocytes. 
In Barnhart et al. (2011), we find the values of viscosity rj ~ 10^ Pa s and active stress ;^co ~ 10^ Pa. The drag 
coefficient can vary over several orders of magnitude depending on the substrate whose physical properties have not 
been specified in Verkhovsky et al. (1999). However, based on the fact thati, in Verkhovsky et al. (1999), the velocity 
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Figure 20: Distribution of myosin (red) and actin (blue) in the static (left) and motile (right) regimes. Insets show the experimental distribu¬ 
tions of actin (cyan) and myosin (red) from Verkhovsky et al. (1999). Picture is taken from http://lcb.epfi.ch/cms/lang/en/pid/71379, courtesy A. 
Verkhovsky. Parameter values: X - 0.0125, = 0.1, = 20 and r = 0.018. 


of the fragment after initiation of motility was approximately 0.08//m s \ we can infer from Fig. 5 of Barnhart et al. 
(2011) that ^ ~ 2 X 10^^ Pa s m“^. From Barnhart et al. (2011) and Luo et al. (2012), we can also obtain the value of 
the diffusion coefficient D ~ 0.25 x 10“^^ m^ s ^ and, from Barnhart et al. (2010); Du et al. (2012); Loosley and Tang 
(2012), we estimate the stiffness of the cortex k ~ 10^ Pa. Finally, directly from Verkhovsky et al. (1999), we infer 
that the characteristic length of the keratocyte fragment is Lq ~ 20 x 10“^ m. Based on these estimates we conclude 
that Z ~ 0.0125, ~ 0.1 and ~ 20. 

In Verkhovsky et al. (1999) (Fig. 5), the initially round fragment with diameter Li = 22 yum, was subjected to an 
applied stress of the order of = 15-20 kPa. The loading was applied after 830 s and lasted for about 80 s. The 
additional surface tractions can be easily incorporated into our model through the boundary condition at the rear of 
the cell: (t(/_( 0, 0 = -[L(0/Lo - q-{t). 

In Fig. 19(a) we present the results of our numerical simulation of the motility initiation experiment of Verkhovsky 
(Verkhovsky et al., 1999). We start with a uniform initial state where motors are distributed homogeneously. We 
chose a generic value of the length L(0) that is slightly different from the value L+ which is unstable in this range 
of parameters. The length first decays towards the value corresponding to the branch as one could expect based 
on Figs. 9(a), 9(b) and 16. This is an unstable state which we found rather robust to selected perturbations. The 
distribution of motors remains non polar with the development of two contractile zones characteristic of the nontrivial 
static regime S ^ • The system then remains in this long living unstable state until we apply an additional one-sided 
force on the boundary breaking the symmetry of the S J state. The destabilized system evolves towards the motile 
state on the branch with both velocity and length well captured by our model. 

We can now compare with experiment the stationary density profiles (for both myosin and actin) generated by the 
model. In the static regime, the flow of actin is absent (v = 0) an the model then predicts uniform distribution of actin 
and myosin. From Fig. 20 (left), we see that this prediction is in agreement with experimental observations given that 
we disregard fluctuations and neglect near-membrane effects. 

From Rubinstein et al. (2009), the turnover time of actin can be estimated to be 30 s. Therefore we obtain in 
non dimensional units that r = 0.018 which leads to the estimate e = 0.015. We recall e accounts for the size 
of polymerization source and sink at the leading and trailing edge, see Section 7 for details. Knowing the value 
of 6, we can reconstruct the mass density distribution p{u) which we show in Fig. 20 (right) together with the motor 
concentration distribution c{u). One can see that outside the boundary layers the model captures the main effect: the 
sweeping of actin towards the de-polymerization zone at the back of the cell by the retrograde flow and its regeneration 
on the polymerization zone at the front of the cell. A more detailed quantitative comparison with experiment requires 
an account of the two (or even three) dimensional nature of the flow. 

Overall, we can conclude that the model reproduces rather well the motility initiation pattern observed in Verkhovsky’s 
experiment. Moreover, the ensuing dynamics is described adequately by the stable motile branch predicted by our 
theory formerly Fig. 19(b). 
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Contractility {V) 


Figure 21: Left: Phase diagram in the parameter plane {P, l/TC) for the system (40) (with no length constraint). The parameter ^ 6 x 10“^ 
is fixed at its experimental value. The solid (red) line indicates the motile bifurcation threshold for the branch Dj" (similar to Fig. 11(b)), while 
the dashed line bounding the metastability domain indicates the location of the turning points on the motile branch in the appropriate analog of 
Fig. 12(b). The dashed line separating static and collapsed configurations indicates the location of the turning point a in Fig. 7. Right: effects of a 
high (top) and low (bottom) concentration saturation thresholds. 


In another experiment by Yam et al. (2007), which we interpret here only qualitatively because of the absence 
of a natural ID representation, motility was induced by injection of calyculin A, known to be a factor increasing 
the activity of myosin motors. The conventional interpretation of this experiment refers to the local variation of 
contractility which disrupts the actin flow and affects the cascade of polymerization and depolymerization (Paluch 
et ah, 2006). Instead, from the perspective of our model it is natural to conjecture that the injection calyculin A 
affects the value of parameter T pushing it beyond the threshold where the static symmetric configuration is stable 
and initiating in this way the polarization instability which in turn leads to motility. 

Notice that in both experiments by Verkhovsky et al. (1999) and by Yam et al. (2007), a fraction of keratocyte 
cells did not move at all after being exposed to the same mechanical or chemical perturbation as the cells that did 
become motile. This can be explained by the fact that the realistic values for P and lay rather close to the boundary 
separating static and motile regimes, see Fig. 19(b). It is then feasible that some cells remain in the symmetric (static) 
regime despite the perturbation. 

It is also feasible that the realistic dependence of active stress on myosin concentration saturates above a certain 
threshold which, as we have seen, can change the nature of the motility initiation bifurcation (Di branch) into a suh- 
critical pitchfork, see discussion in Section 5. This opens a finite range of metastability where both the homogeneous 
static state and the inhomogeneous motile state are locally stable. The implied non-uniqueness may be an alternative 
explanation of the simultaneous presence of motile and non motile cells despite apparently equal levels of contractility. 

To exemplify this last claim, we show in Fig. 21 the effect of switching to threshold type dependence of contractile 
stress on the concentration of motors, see (40). Notice, that we have dropped in Fig. 21 the assumption that the length 
of the moving segment is fixed. A comparison of Fig. 21 with Fig. 1 1(b) shows that the saturation of contractile stress 
introduces a finite zone of metastability between static and motile configurations: in this zone finite perturbations are 
required to switch from static to a motile regime. This prediction was very recently confirmed in vivo by Barnhart et al. 
(2015) and the metastability domain as in Fig. 21 was mapped experimentally. We also observe that for sufficiently 
large values of the active stress saturation threshold r, our model associates metastability with both, motility initiation 
and motility arrest. On the arrest side (Lancaster and Baum, 2014), this prediction can be finked to the metastability 
of cell division (Turlier et al., 2014) which to our knowledge has not been yet experimentally documented. 
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9. Conclusions 


We studied a prototypical model of a crawling segment of an active gel showing the possibility of spontaneous 
polarization and steady self propulsion in the conditions when contraction is the only active process. Our model, 
which focuses entirely on ‘pullers’, complements the existing theories of polarization and motility that place the 
main emphasis on ‘pushers’ and link motility initiation with active treadmilling and protrusion. Mathematically, the 
proposed model reduces to a dynamical system of Keller-Segel type, however, in contrast to its chemotaxic analog, 
the nonlocality in this model is due to mechanical rather than chemical feedback. If compared with previous studies 
of Keller-Segel type problems, our setting is complicated by the presence of free boundaries equipped with Stefan 
type boundary conditions. 

As we argue, the motor proteins with sufficient contractility induce internal stress which can overcome the hydro- 
dynamic resistance and induce flow. The flow produces a drift of motors in the direction of the regions where they 
concentrate and such autocatalytic ampliflcation is the mechanism of the positive feedback in our model. The ensuing 
runaway is countered by diffusion of motors which penalizes creation of concentration gradients and thus plays the 
role of a negative feedback. When a critical contractility of motors is reached, the homogeneous distribution of motors 
becomes unstable. The contraction asymmetry then induces a flow of actin filaments towards the trailing edge thus 
producing frictional forces which propel the cell forward. The rebuilding of the balance between drift and diffusion 
leads to the formation of a pattern. Among various admissible patterns, whose number increases with contractility, 
the stable ones localize motors at the trailing edge as observed in experiments. 

The proposed model provides an alternative qualitative explanation of the experiments of Verkhovsky et al. (1999) 
and Yam et al. (2007) that have been previously interpreted in terms of active polymerization inducing the growth of 
actin network (Blanch-Mercader and Casademunt, 2013). Most strikingly, the predictions of our model are also in 
quantitative agreement with experimental data, which is rather remarkable in view of a schematic nature of the model 
and the absence of fitting parameters. 

In addition, the model captures a durotactic effect since the directional motion cannot be initiated if friction with 
the substrate is larger than a threshold value. We show that below this threshold, motile regimes exist in a finite range 
of contractility. This means that if the cell is already in motion, it can recover the symmetric (static) configuration 
either by lowering or by increasing the amount of operating motors. The predicted possibility of cell arrest under the 
increased contractility should be investigated in focused experiments. 

We have also shown that when the contractility depends on the motor concentration nonlinearly, one can have a 
metastability range where both static and motile regimes are stable and can coexist. In this range of parameters a 
mechanical perturbation may be used to switch back and forth between static and dynamic regimes and reproducing 
such behavior in vivo presents an interesting challenge. This prediction of the model is particularly important in the 
context of collective cell motility (in tissues) where contact interactions are able to either initiate or terminate the 
motion (Abercrombie, 1967; Heckman, 2009; Trepat et al., 2009; Vedula et al., 2012). 

Despite the overall success of the proposed model, it leaves several important questions unanswered. Thus, our 
focus on a one dimensional representation (projection) of the motility process obscured the detailed description of the 
reverse flow of actin monomers which we have replaced with an opaque jump process. Similarly, our desire to max¬ 
imally limit the number of allowed activity mechanisms, forced us to assume that polymerization of actin monomers 
and their transport are equilibrium processes. The assumption of infinite compressibility of the cytoskeleton, which 
is behind the decoupling of the mass transport from the momentum balance, also remains highly questionable in 
the light of recent advances in the understanding of cytoskeletal constitutive response (Broedersz and MacKintosh, 
2014; Pritchard et al., 2014). Finally, our schematic depiction of focal adhesions as passive frictional pads needs to 
be corrected by the account of the ATP driven integrin activity and the mechanical feedback from the binders to the 
cytoskeleton (Schwarz and Safran, 2013). These and other simplifying assumptions would have to be reconsidered 
in a richer setting with realistic flow geometry which will also open a way towards more adequate description of 
membrane (cortex) elasticity and will allow one to account for the polar nature of the gel (Marchetti et al., 2013). 

Ultimately, the answer to the question whether the proposed simplified mechanism provides the fundamental 
explanation of motility initiation in keratocytes will depend on the extent to which the inclusion of all other related 
factors affects the main conclusions of the paper. A more thorough analysis will also open the way towards deeper 
understanding of the remarkable efficiency of the proposed mechanism of self-propulsion delivering almost optimal 
performance at a minimal metabolic cost (Recho et al., 2014). 
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Appendix A. Analysis of nontrivial static solutions. 

Solutions of boundary value problem (29) correspond to closed trajectories on the phase plane (s, s') passing 
through the origin (s = 0, s' = 0) and different types of such trajectories are illustrated in Fig. 22. 



Figure 22: Phase diagram for the static solutions in the parameter space (A, B). A + B = 0 line is the trivial (homogeneous) solution. In the bottom 
corner we show the blow up of the same diagram. 


Depending on the position of a point in the parametric plane A, B, one can identify five different types of behavior: 


1. If A + .S = 0, then equation TE(r) = 0 has one double root at r = 0 and one single root (negative or positive) at 
r = s- (Case 3 on Fig. 22). The only solution is then the trivial one s(u) = 0 and L = L+. 

2. If A + .S < 0, then Eq. W(r) = 0 has three roots: r = 0, r = 5'_ < 0 and r = 5'+ > 0. This case corresponds 
to static branches labeled in Fig. 22, Fig. 7, Fig. 9(a) and Fig. 9(b) by numbers without a prime (Case 1 on 
Fig. 22). In this domain we find nontrivial static solutions with 0 < s(u) < s+. Different solutions correspond 
to different number (m) of sign changes for the function s'(u) and different values of L = 2m W{cr)~^dcr. 


3. If 


A + 5>0 

1 - Va2 -2B+\ < 5e-VA2-2fi+i+A+i (A.l) 

A < -1 


then Eq. W{r) = 0 has three roots: r = 0 and r = 5’_ < 0 and r = s+ < f) with, 5’+ > 5’_. This case 
corresponds to non motile branches labeled in Fig. 22, Fig. 7, Fig. 9(a) and Fig. 9(b) by numbers with a prime ' 
(Case 2 on Fig. 22). In this domain, we find nontrivial static solutions with S- < s{u) < 0. Again, different 
solutions correspond to different number of sign changes for the function s'{u) and different values of L = 
2m W(cr)~^dcr. 
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4. If 


A + 5>0 

^ 1 - Va 2 -25+1 < 5^-VA2-2fi+l+A+l (A.2) 

A > -1 

then Eq. iy(r) = 0 has three roots: r = 0 and r = i-. > 0 and r = s+ > 0 with, 5-+ > 5'_ and there are no static 
solutions since there are no closed paths in the phase plane passing through the point (0,0). 

5. If 1 - Va^ -25+1 > Be~ Va2-25+i+a+i^ equation W(r) = 0 has only one non degenerate root Situ = 0. In 
this case there are no static solutions since again there are no closed paths in the phase plane. 

Notice also that for the solutions described above the map between the two parameterizations (A, B) and (7C, P) is 
explicit 

= A/(2m £* W(cr)-^^^d(r - 1 ) 

TCP = 2m — A)W{<T)~^^^d<T. 

Notice also all nontrivial static solutions bifurcate from the trivial branches in the sense that there are no detached 
solutions. Indeed, if a solution were detached, it would not pass through the origin (trivial solution) in the space (s, s'). 
But that would mean it cannot satisfies boundary conditions. 


Appendix B. Analysis of the characteristic equation (39). 

Equation (39) has two families of solutions depending on whether oj is real or pure imaginary. In the first case, we 
denote ojc = |u;| > 0 whereas ojc = -\oj\ <0 in the second case. 

2[cosh(6t;c) - 1] + (ZojI/L?' - \)cOc sinh(u;c) = 0 if > 0, tR 1 i 

2[cos(ojc) - 1] + (ZcollL^ + l)oJc siniojc) = 0 if u;? < 0. ^ ^ 

Equations (B.l)i and (B.l )2 need to be analyzed separately: 

1. When > 0, equation (B.l)i has a unique solution only if L^/Z > 12. It is given by the implicit formula 

2tanh(y) = (1 - 


The unstable eigen-vector can be written as 


' 6L " 


^ 0 

6V 

= 

1 

, ds(u) ^ 


1/2)0;,] (2 m l)sinh(o;,/2)} J 


Since SV ^ 0, this bifurcation leads to a motile configuration which we denote Di . In Eig. 23 the eigen¬ 
functions associated with the sub-branch Z)| bifutrcating from the trivial solution L+ are illustrated for Z = 
0.01. As parameter Z increases the exponential viscous boundary layers thicken. They fully disappear at 
Z = T^/12 where the ‘hyperbolic’ eigen-vectors become ‘trigonometric’. 

2. When < 0, equation (B.l )2 has two sub-families of solutions: 

(a) The first family can be written explicitly : cl>c = -Imn with m > 1. The unstable eigen-vector has the form 


^ 6L ' 


6V 

= 

, Ss(u) , 

V 


1 

0 

{2L-1){ZojI+P) Y / X n 

LKt-h 1 ] 




/ 


Since dV = 0, the bifurcated branch describes the nontrivial static solutions Sm studied in Section 4. In 
Eig. 23 the eigen-functions associated with the sub-branch S'^ originating at the trivial solution L+ are 
illustrated for = 0.15. 
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Figure 23: Solution branches of the characteristic equation (B.l) as functions of Z. for the trivial static solution L+ (!P = 0.01). From (36), the 
locus of the bifurcation points are recovered and shown of Fig. 8 . We refer to Fig. 9(a) for the label of bifurcation points. We represent in inserts 
the eigenfunctions 6s related to ^J, £> 3 , 3 for = 0.15 and the eigenfunction 6s related to Dj" for Z = 0.01. The eigenfunctions are 
normalized to 1 ; solid and dashed lines correspond to the two possible directions of the pitchfork bifurcation. 


(b) The second family consists of a countable set of negative roots of equation (B.l )2 given implicitly by, 

2tan(y) = (B.2) 

The unstable eigen-vector is 


' (5L ' 

f 

6V 

= 

. Ss{u) , 

V 


V Z<^c COS(6t>c/2) 


0 

1 

sin[(w - l/2)6t;c] - 2{u - 1/2) sin(6G)c/2)} 


It corresponds to motile branches because 6V 0. We denote this family by D^. In Fig. 23 the eigen-functions 
associated with a subbranch originating at trivial solutions L+ are illustrated for X = 0.01. 


Appendix C. Normal forms. 

Normal form in 7C. In terms of the normalized stress variable r = s/L the original nonlinear problem can be written 
as 


Mr(u)-Vu) 


-Zr'\u) -f L^r(u) = TCP— 


with r(0) = r(l) = 0, r'(0) = r'(l) = V. 


JT gL(r(u)-Vu)^^ 

Assume that 6 is a small parameter and expand the solution of (C.l) around a bifurcation point up to third order 


(C.l) 


2 3 

2t//o , ^3t 


' 2^0 . .3? 


r = 0 -f 6r -f 6^/2 -f 6^/6 -f 0 ( 6 "), V = 0eVjl -f e^VI6 + o{e% L = L + eL + e^L/l -f 6"L/6 -f o{e^). 
Assume that the bifurcation parameter TC and therefore 


2 3 


+ 67C + e^TCjl -f e^TCje -f o{e% 
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1 2 
Figure 24: (a) Values of TC for the S and branches as functions of parameter ^ for = 0.001. (b) Values of TC for the and D~ branches 

2 

as functions of parameter !P for Z = 0.001. The point where = 0 along the L_ static branch indicates a nature change of the motile pitchfork 
bifurcation from supercritical to subcritical pitchfork. The parameter dependence of such a point is represented as a function ofP = in the 

inset. 


0 

where is the bifurcation point. These expressions are then inserted into equation (C.l). Separating different orders 
of 6 we obtain three differential equations 

0(1), L(r,i,y) = 0, (C.2) 

0(2), L(r, L,V) = ^ Po(^ L, V) + '^ Pi (r, L, V), (C.3) 

0(3), L(r, L, y) = 7C Qo(s, L, V, r, L,V)+'K Qi(r, L, V, r, L,V) + ‘K Q 2 (r, L, V, I, L, V), (C.4) 


where L is the linear operator already introduced in the stability analysis, see (34): 
L(r, L,V) := r"{u) - (j/r(u) + 




Z (L-l)L^ 

and Pq, Pi, Qo. Qi and Q 2 are known non linear operators. The boundary conditions remain the same at all orders i: 


r(0) = r(l) = 0, r'(0) = r'(l) = V 

0 

In the leading order, we obtain the results already reported in Section 5 including the eigenvalue and the 
1 1 1 

eigenfunction r(u), L, V. To have a nontrivial solution in the next order, the right-hand side of equation (C.3) must be 
orthogonal to the kernel of the dual of L (for the scalar product). In the (Ci, C 2 , SL, SV) space, see Section 5, this 
means orthogonality to the kernel of the transpose of matrix (38). The resulting linear scalar equation determines the 

value of When this value vanishes, higher orders must be considered in a similar way. We summarize below the 
main results obtained by implementing this procedure. 

1. Static branches result from transcritical bifurcations. For the m* branch we have 

‘k = {C- Am^n^Z)l[C{L - 1 )] 

1 

2. Motile branches all correspond to pitchfork bifurcations with = 0. They can be either subcritical or super- 
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critical depending on the sign of 

2 

7C = ((2L - (?>a? + 770) + + 1415(1 - 2L)a? + 4620(1 - 2L)) + + 610;^+ 

374) - 31^"^ - 1340a;2 - 7480) + 2L^cj'^Z^ {al(6cjJ^ + - 3150) - 9cJ^ + 50^^ + 7380) + [-2L(165cJ^+ 

65()2a? + 23160) + l95co^ + 6574^^^ ^ 22440) + 6tco^Z^ ((61L - 34)0;"^ + (2622L - 1339)^^ + 7072(2L - 1)) + 

L^oj^^Z^ (3(31 - 60L)aj'^ + 4264(1 - 2L)a? - 28224(2L - 1)) + 2{2L - l)a^^ {9aj^ + A12u? + 3456) Z^) 

(l44(L - 1)(2L - l)a;^Z^ [l^ - a?Z) - 2L^ + 6)Z + aj^Z^)y^ (C.5) 


2 2 

In Fig. 24(b) we illustrate the function 7C(^) for the first motile branches D~ and Z)[. As 7C > 0 for all values of 
the activity parameter P, the motile branch Z)| always bifurcates from the static branch in a supercritical (pitchfork) 
manner. In contrast, the motile branch D~ can bifurcate either supercritically or a subcritically depending on the value 

2 

of P. When P is larger than a threshold value Ps, the coefficient 7C changes sign and becomes negative indicating a 
subcritical character of the bifurcation on the L_ static branch. 


Normal form in P. We now consider P as the bifurcational parameter. The derivation of the normal form in this case 
is more complex because the homogeneous static solution LiP) is a multivalued function of P (see Fig. 6). One can 
circumvent the difficulty by introducing a new variable 


/ = L(L - 1) + P, 

Then the boundary value problem (C.l) takes the form 


(C.6) 


^L{r{u)-Vu) 


-Zr"(u) + L^r(u) = %[J - L(L - 1)] — 




, with r(0) = r(l) = 0, r'(0) = r'(l) = V. 


whose trivial solution is (/, V, r) = (0,0,0). In this formulation /, V and r(u) are unknowns while the length L is the 
bifurcation parameter. The regular expansions near the homogeneous state give 

r = er + e4/2 + e^r/6 + o(e^), V = eV + + e^y/6 + o(e^), J = eJ + jl + ed/6 + o(e^), 


L = L + eL + efl + ed/6 + o(e^). 

Distinguishing the static and motile branches as before, we obtain the following results: 


0 


1. Static branches are all found to be transcritical bifurcation. For the m* branch, we have P = P eP + o{e) 
where 

0 0 0 1 0 1 
P = -L{L -\),P=\-{2L - 1)L. 


and L is a solution of the cubic equation 


\2 ° 


-'K{Lf{L - 1) = + {Lf. 


In this equation only two roots corresponding to 5^ (the smaller) and to 5^* (the larger) are in the range [0,1]. 

1 

In Fig. 25(a), we illustrate the behavior of of the function PifK) for the branches S ^ and S 
2. Motile branches result from pitchfork bifurcations that can be either supercritical or subcritical with 




p = p + e^Pjl + o{e^). 
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Figure 25: (a) Values of V for first static branch as a function of parameter TC for a fixed = 1. (b) Values of P for the first motile branch {X = 1). 


The coefficients in this expansion can be written inn the form 

0 0 0 2 2 0 2 
r=-L{L-\) and r = J-{2L-\)L, 

where 

2 _ 60LZ - TCCL - l)(i)3(7C(L - 1) - 4) 

^ ~ 0 • 

2A%{%{L - 1) + 1)2 

0 

The length L can be found from the system of equations 

0.0 . 0 ^ 

- 1) = + (L)2 

tanh(w/2) = (w/2)(l - Zu?l{Lf) 

Again, two roots are in the interval [0,1]: the smaller one belongs to the branch and the larger one to 

2 

the branch In Fig. 25(b), we illustrate the function for m = 1. The bifurcation from the static 

2 

homogeneous solution with longer length is always supercritical as V(JT^ > 0. Instead, the bifurcation from the 

2 2 

Static homogeneous solution with smaller length can change from subcritical < 0) to supercritical (^ > 0). 


Appendix D. Normal form for the stiff limit. 

In the study of (41) we closely follow the procedure developed in Section 5. In essence, the results are exactly 
the same for fixed L = 1 and the product replaced by A with only one homogeneous state (s(y) = 0, V = 0 and 
So = 0(r)/r) and where Ssq replace SL. As a result, there is an infinite sequence of bifurcations branching from the 
now unique homogeneous state. We shall only focus on the stable attractor of the problem, namely, the homogeneous 
solution before the Di bifurcation and the first motile branch after. 

The critical value of the bifurcation parameter Ac corresponding to the case when a homogeneous static solution 
becomes linearly unstable is given by the formula Ac = (1 + r)^(l - ZoJ - c^)- where ojc is a root of the equation 
tanh(6L)c/2) = ojc(l - X(j^^c)I^ with the smallest absolute value. We then proceed to the next order developing a regular 
expansions close to the bifurcation point 


s = 0 6S + 6^i/6 + o{e^) 

y = 0 + 6y + 62y/2 + 6^y/6 + o{e^) 

12 3 

So = 0(r)/r + eso + €^sol2 + €^sol6 + o(€^), 
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Figure 26: Parameter A characterizing the structure of the static to motile bifurcation in the case of non linear contraction law. 


Similar expansion for the bifurcation has the form 


: Ac + 6/1 + 6 ^/ 1/2 + ^Aj^ + ( 9 ( 6 ^). 


For the first motile branch one finds that A = 0, indicating a pitchfork bifurcation. 

Below we show that this bifurcation can change from supercritical to subscritical depending on the value of the 

1 2 

dimensionless parameter r. Assuming without loss of generality that V = V = 1, we obtain 

(o/Z-l){Ar^ + Br + C) 


2 

A = 


144oj^Z^ {oj^Z^ - 2 (u;2 + 6) Z + 1) 


where 


A = - l23aj^°Z‘^ + 60^(35 - 164Z)Z^ + 2w®Z^(1073Z - 84) + 6w4z(l440Z^ - 155Z + s) 

+3o/ (l320Z^ - 430Z + l) - 9240Z + 770 

B = -2(2W^Z^ - 87w'0z4 + A(J’{39 - 173Z)Z^ + 2w®Z^(707Z - 66) + 3oj'^Z{l440Z^ - 210Z + 13) 

+0? (2280Z^ - 1150Z + 3) - 9240Z + 770) 

C = -6oj^^Z^ + 2loj^^Z^ + 6w*Z^(28Z - 5) + 2w®(12 - 79Z)Z^ + 6oj^Z(S5Z - 2) + 3oj^ (l320Z^ - 430Z + l) 

-9240Z + 770 


these expressions show that there exists a critical value of the parameter r such that the bifurcation is supercritical 

2 

(i.e. /I < 0) for r < This regime corresponds to a state where contraction is proportional to concentration of motors. 

2 

For r > Vc, the picthfork bifurcation is subcritical (i.e. /I > 0) and the regime is characterized by a contraction which 

2 

saturates into the plateau. We plot on Fig. 26 the value of /I as a function of r for a fixed Z = 1 and in inset the value 
rdZ). When Z ^ 0, ^ 2 and when Z ^ oo, ^ (7 + V^)/10. 
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